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DEFENSE  METIEORCEiOGICftL  SATELLITE  PROGRAM  PROCESSING 


1.0  INTRODUCTION 

1.1  Background 

BEDFORD  RESEARCH  Associates'  work  on  the  "J4",  "J*",  "IE",  and  "M" 

experiments  of  the  Defense  Meteorological  Satellite  Program  (DMSP)  which 
consists  of  processing  raw  data  tapes,  and  creating/maintaining 
applications  programs.  The  data  from  the  F2  through  F7  satellites  is  put 
onto  DM  tapes  at  Offut  Air  Force  Base,  Nebraska.  These  tapes  are  shipped 
to  the  Air  Force  Geophysics  Laboratory  (APGL)  for  processing  by  BEDFORD 
RESEARCH  Associates  and  the  Space  Physics  Division.  Each  experiment's  data 
is  first  processed  by  a  program  designed  to  create  a  database  (this  is 
Phase  1) .  Phase  2  consists  of  creating  interactive  programs  to  use  this 
data.  BEDFORD  RESEARCH  Associates  designs  these  programs  and  maintains  the 
"J4"  and  "J*"  databases. 

2.0  PROJECT  DESCRIPTION 

Each  experiment  has  a  database  interactive  package  and  an  applications 
interactive  package.  Each  experiment's  database  package  has  the  following 
essential  options:  1)  an  option  to  merge  ten  DM  tapes  onto  one  tape,  and  2) 
an  option  to  take  up  to  five  of  the  latter  and  reduce  the  data  to  a 
standard  format.  The  latter  process  is  called  a  production  run.  Other 
options  are  also  present  in  the  database  packages,  depending  on  the 
experiment.  The  applications  programs  for  each  experiment  are  similar. 
The  following  describes  the  interactive  packages  in  more  detail. 

2.1  " J 4"  Database  Interactive  Package 

The  procedure  which  runs  the  "J4"  database  interactive  package  is  J4INTCR. 
When  this  command  is  issued,  the  user  will  be  prompted  with  a  series  of 
questions  that  enable  him  to  create  one  of  the  following  batch  jobs: 

1)  Merge  ten  DM  tapes  to  one  tape 

2)  Create  a  database  fcy  using  up  to  five  of  option  l's  tapes 

3)  Create  a  Computer  Center  (CC)  tape  of  data  from  up  to  five  of  option 
2's  tapes 
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Backups  of  the  CC  tapes  are  also  made  by  BEDFORD  RESEARCH  Associates 
personnel  and  kept  at  the  base  office. 

The  tapes  from  option  1  are  given  a  magnetic  label  of  HAX  followed  by  three 
lumbers.  The  paper  label  is  BR  followed  by  three  numbers.  All  tapes, 
except  the  CC  tapes,  are  supplied  by  BEDFORD  RESEARCH  Associates.  These 
tapes  are  reduced  to  standard  format  by  half-months  for  each  satellite  (F6 
and  F7)  by  using  option  2. 

Option  2,  a  production  run,  prompts  the  user  for  the  satellite  desired.  It 
also  prompts  the  user  for  the  beginning  or  end  of  a  month.  The  range  for 
the  beginning  of  the  month  is  day  1  through  day  15;  the  range  for  the  end 
of  the  month  is  day  16  through  day  31.  The  tape  is  given  a  magnetic  label 
of  br  as  well  as  a  paper  label  of  BR.  An  example  of  a  multifile  label  on  a 
production  run  tape  is  B1AFR85F6  (beginning  of  April  [1-15] ,  1985, 
satellite  F6) . 

These  production  run  tapes  are  then  copied  to  a  CC  tape  by  option  3.  These 
CC  tapes  contain  up  to  five  half-months  of  data.  For  example,  the  aid  of 
February,  1985  through  the  end  of  April,  1985  for  satellite  F7  might  be  on 
one  CC  tape.  The  data  for  F6  and  F7  production  runs  are  kept  separate. 
These  CC  tapes  are  used  by  the  Space  Physics  Division  for  the  applications 
interactive  par’  age.  They  are  also  backed  up  by  BEDFORD  RESEARCH  personnel 
in  case  a  CC  tape  is  destroyed.  For  a  brief  outline  of  the  components  of 
the  J4  database  interactive  package,  see  Appendix  A. 

2.2  "J*"  Database  Interactive  Package 

The  "J*"  database  interactive  package  is  similar  to  the  one  for  "J4".  Ten 
DM  tapes  are  merged  onto  one  tape  with  a  magnetic  label  of  STIR  followed  by 
three  numbers.  The  paper  label  is  a  BR  label. 

These  tapes  are  then  reduced  to  standard  format  and  put  onto  an  OS  tape. 
The  OS  number  is  supplied  by  the  computer  center.  Unlike  "J4",  a  whole 
month  of  data  is  put  on  an  OS  tape.  Also,  wily  the  F7  satellite  is 
involved. 
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These  OS  tapes  are  further  reduced  by  packing  the  data.  The  data  on  one  OS 
tape  is  packed  and  put  onto  another  OS  tape.  These  OS  tapes  are  used  by 
the  Space  Physics  Division  in  the  applications  programs.  These  tapes  are 
also  backed  up  by  BEDFORD  RESEARCH  Associates.  For  an  outline  of  the  "J*" 
database  package,  see  Appendix  B. 


WW.TOT 


2.3  "IE"  Database  Interactive  Package 

The  "IE"  database  interactive  package,  used  by  the  Space  Physics  Division, 
has  five  options.  The  first  option  merges  the  data  and  the  second  reduces 
it  onto  DMM  tapes  by  half-months.  Option  3,  however,  used  the  production 
runs  to  create  four  files  of  reduced  data:  Ion  Sweep,  Electron  Sweep,  Ion 
Density,  and  Electron  Density.  Option  4  allows  one  to  copy  the  information 
from  Option  2  or  3  onto  a  tape  for  use  in  the  applications  package.  Option 
5  copies  a  tape  for  backup  purposes.  See  Appendix  C  for  an  outline  of  the 
components  of  the  "IE"  database  package. 

2.4  "M"  Database  Interactive  Package 

The  "M"  database  package  as  three  options.  The  first  merges  ten  raw  data 
tapes  onto  one  tape.  The  second  reduces  these  to  one  DMM  tape.  This  tape 
is  then  copied  to  a  CC  tape  and  backed  up.  This  package  is  very  similar  to 
the  "J4"  package.  However,  two  versions  of  the  "M"  package  are  used  by  the 
Space  Physics  Division  —  one  for  data  prior  to  January  13,  1985,  and  one 
for  data  for  Jan.  13  -  Jan.  31.  The  data  format  changed  to  ASCII  and 
subsequently  the  packages  were  changed.  See  Appendix  D  for  an  outline  of 
the  two  versions. 

2 . 5  "J 4". Applications  Interactive  Package 

There  are  four  options  in  this  package:  list,  flux,  spec,  and  file.  List 
creates  a  listing  of  the  requested  data.  Flux  creates  flux  plots  of  ion 
and  electron  data.  Spec  creates  spectra  plots,  and  file  creates  a  file  of 
distr itwtion,  differential  flux,  or  counts  data.  All  options  can  be  run 
for  whole  day  intervals  or  specific  time  periods  of  less  than  a  day.  Also, 
if  one  chooses  the  flux  option,  one  can  obtain  a  printout  of  the  values  of 
JTOT,  JETOT,  and  EAVE. 

These  options  use  the  CC  tapes  created  by  the  database  packages.  The 
procedure  file  creates  a  batch  job,  which  is  a  local  file.  This  file  can 
then  be  routed  by  the  user.  See  Appendix  E  for  an  outline  of  the  package. 
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2.6  "j2*"_Appiications  interactive  _ package 

There  are  seven  options  in  the  "J*"  package.  Each  can  be  run  for  intervals 
in  whole  days  or  for  smaller  time  periods.  The  options  include  the 
following:  1)  list,  2)  print  sums  of  counts,  doses  or  fluxes,  3)  file,  4) 
survey  plots,  5)  output  non-zero  star  counts,  6}  polar  pass  averages,  and 
7)  polar  plots.  See  Appendix  F  for  an  outline  of  the  package's  components. 

2.7  " i£ll.j\ppiicatioris_interftctivft.  Package 

This  is  very  similar  to  the  "J4"  package.  It  has  three  options:  list, 
flux,  and  file.  There  is  no  spectra  optionn.  See  Appendix  G  for  an 
outline  of  the  components  of  the  package. 

2.8  "M"  Applications  Interactive  Package 

This  is  also  very  similar  to  the  J4  package.  However,  there  are  five 
options:  list,  flux,  spec,  file,  and  magf.  The  latter  generates  plots  of 
the  magnetic  files.  See  Appendix  H  for  an  outline  of  the  components  of  the 
package. 

2.9  SunBary_pfJEackages 

Basically,  there  are  four  data  base  packages  and  four  applications 
packages.  BEDFORD  RESEARCH  Associates  maintains  two  of  the  data  bases. 
The  applications  programs  are  used  by  the  Space  Physics  Division  of  AEGL. 
They  also  maintain  the  "M"  and  "IE"  databases.  All  the  packages  are  for 
the  F6  and  F7  satellites,  except  for  the  "J**  package  which  is  for  the  F7 
only. 

3.0  PROCESSING  ACCOMPLISHED 

Data  processing  has  been  done  on  the  following: 

IE  -  from  November  1983 

J*  -  from  November  1983 

J4  -  from  December  1982 

M  -  from  November  1983 


4.0  AEDITIONAL  EFFORTS 


It  had  also  been  requested  that  a  file  option  for  "J4"  be  created  to 
contain  the  ephemeris  and  data  for  the  electrons,  channels  1-13,  for 
satellite  F7.  Also,  the  values  JETOT,  JTOT,  EAVE,  £  and  £  0  have  been 
calculated,  printed,  and  put  in  a  file.  The  following  are  the  formulas  for 
the  last  three  values: 

JETOT 

EAVE  =  -  »  EQ 

JTOT 

1/2  20EO 

Sp  =  I  - ,  where  IE  =  2  *  1.602  x  10~10JETOT 

E  (4+Eq2) 


5/8 

“  H  =  Eo  -  p 
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APPENDIX  A 


HLIBCOD 


J4INTCR: 

INTERAC: 

LIBJ4CR: 

HLIBCOD: 


"J4"  DATABASE  INTERACTIVE  PACKAGE 


J4INTCR 


INTERAC 


LIBJ4CR 


Procedure  file  which  runs  the  package 

FORTRAN  program  which  creates  the  batch  file 

Library  of  FORTRAN  code  which  merges,  reduces,  and 
copies  the  data. 

FORTRAN  code  which  merges,  reduces,  and  oopies  the 
data. 
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APPENDIX  B 


JSIRCR 


JSINTCR: 

JSINTER: 

JSLIBCR: 

JSTRCR: 


"J*"  DATABASE  INTERACTIVE  PACKAGE 


JSINTCR 


JSINTER 


JSLIBCR 


Procedure  which  runs  the  package 
FORTRAN  code  which  creates  the  batch  job 
Library  code  of  JSTRCR 

FORTRAN  code  which  merges,  reduces,  packs  snd 
copies  the  data 
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APPENDIX  C 


MLIBCOD 


IEINTCR: 

MINTERA: 

IELIBCR: 


"IE"  DATABASE  INTERACTIVE  PROGRAM 


IEINTCR 


MINTERA 


IELIBCR 


Procedure  to  run  the  package 
FORTRAN  code  which  creates  a  batch  job 
Library  of  code  of  MLIBCOD 

FORTRAN  code  which  merges,  reduces,  and  copies  the 


MLIBCOD: 


APPENDIX  D 


"M"  DATABASE  INTERACTIVE  PACKAGES 


Befoie_ASCXI 

MDICXJLD 

i 

OLEMIC 

MDLIBCO  — ►  MDLIBCX* 

*  uses  magfldd 
**  uses  IGRF80 


ASCII 

(©INTER 

l 

(©INTER 

MDL3BCR  —  MDLIBC1** 


MDICOLD,  (©INTER: 
CJLDMIC,  (©INTER: 
(©LIBCX,  MDL3BC1: 
(©LIBCO,  MDLIBCR: 


Procedures  which  rw  the  packages 
FORTRAN  programs  which  create  batch  jobs 
Libraries  of  FORTRAN  code 

FORTRAN  codes  which  merge,  reduce,  and  copy  data 
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APPENDIX  E 


"J4"  APPLICATIONS  INTERACTIVE  PACKAGE 


J4INTAP 


t 

HARDYTA - - - -  HAFLE 

_ l_ 


HAPPCODt  - - -  LIBJ4A1  LIBJ4A2  - -  HBCODE 


J4INTAP: 

HAFLE: 

HARDYTA: 

LIBJ4A1 : 

LIBJ4A2: 

HAPPCOD: 

HBCODE: 


Procedure  which  runs  the  package 

FORTRAN  code  which  creates  the  batch  job 

List  of  data  tapes  available  for  use 

Library  of  code  for  list,  flux,  and  spec  options 

Library  of  code  for  file  option 

List,  flux,  and  spec  FORTRAN  code 

File  FORTRAN  code 


APPENDIX  F 


"J*"  APPLICATIONS  INTERACTIVE  PACKAGE 


JSINTAP 


KRJSTAP  - -  J  STRAP 


JSFORAP  - - -  jaiRAP 


JSINTAP:  Procedure  to  run  the  package 

J  STRAP:  FORTRAN  code  to  create  batch  job 

KHJSTAP:  List  of  tapes  available 

JSLIBAP:  Library  of  code 

Code  to  produce  the  data  for  plotting,  etc 


JSFORAP 


APPENDIX  G 


"IE"  APPLICATIONS  INTERACTIVE  PACKAGE 


MULLENS 


MAPFCOD 


IE INTAP 


MAPLE 


IELIBAP 


IE INTAP:  Procedure  which  runs  package 

MAPLE:  FORTRAN  code  which  creates  batch  job 

MULLENS:  List  of  tapes  available 

IELIBAP:  Library  of  code 

Code  which  processes  data  for  an  option 


MAPFCOD 


APPENDIX  H 


MDHAPFC 


MDINTAP: 

MDHAPLE: 

MELIBAP: 


"M"  APPLICATIONS  INTERACTIVE  PROGRAM 


MDINTAP 


MOHAHjE 


MDLIBAP 


Procedure  which  runs  package 

FORTRAN  code  which  creates  the  batch  job 

Library  of  code 

Code  to  process  data  for  options 


MDHAPFC: 


Cn2  CALCULATIONS  FROM  RAWINSCNDE  DATA  USING  THE  VANZANDT  AND  KAIMAL  MODELS 

1.0  INTRODUCTION 

1.1  Overview 

C2,  the  refractivity  structure  constant,  is  the  most  important  parameter 
in  determining  the  effects  of  clear-air  turbulence  on  optical  beam  propag¬ 
ation.  Various  optical  properties,  such  as  the  transverse  coherence  length 
and  the  isoplanatic  angle,  are  determined  from  altitude  weighted  integrals 
of  Cn2  through  the  atmosphere  down  to  the  observation  point.  1,2,4  The 
dominant  contribution  to  these  integrals  will  usually  be  from  the  atmos¬ 
pheric  boundary  layer  but  the  contribution  from  the  free  atmosphere,  espe¬ 
cially  in  the  presence  of  a  strong  jet  stream  is  not  negligible.  Also,  the 
existence  of  a  strong  temperature  inversion  layer,  marking  the  transition 
from  the  boundary  layer  to  the  free  atmosphere,  can  contribute  signifi¬ 
cantly  to  the  integrals. 

1.2  Background 

Various  methods  exist  for  obtaining  altitude  profiles  for  Cn2,  such  as 
thermosonde  radar  and  stellar  scintillometers.  To  develop  seasonal  pro¬ 
files  at  various  locations  throughout  North  America  which  could  be  used  in 
best/worst  case  studies,  for  example,  it  is  desirable  to  develop  methods  of 
calculating  Cn2  from  raw insonde  measurements  of  temperature,  pressure  and 
wind  velocities.  The  main  advantage  is  the  extremely  low  cost  per  raw  in¬ 
sonde  launch  as  compared  to  other  more  direct  but  expensive  techniques  used 
to  determine  Cn2.  Models  developed  using  the  larger  raw insonde  data  base 
can  be  compared  to  measurements  from  these  other  techniques  to  determine 
reliability.  In  the  application  of  raw  insonde  data  to  obtain  Cn2  esti¬ 
mates,  we  have  made  use  of  two  models,  one  for  the  free  atmosphere  and  the 
other  for  the  boundary  layer  to  calculate  Cn2  from  0.1  km  up  to  4  km. 

2.0  DATA 

For  the  model  development,  we  have  used  measurements  from  standard 


raw  insondes  that  were  modified  to  telemeter  temperature  data  at  a  hi^ier 
rate  than  is  normally  obtained  and  required  for  weather  reporting  purposes. 
The  higher  resolution  is  required  here  in  order  to  determine  accurately  the 
first  temperature  inversion  layer  above  the  ground  for  purposes  explained 
further  in  discussions  of  the  models  used.  These  data  were  obtained  by 
personnel  at  the  Pacific  Missile  Test  Center,  California. 

The  measurements  were  made  at  3  sites  within  75  miles  of  one  another  in  the 
proximity  of  Los  Angeles,  California.  A  total  of  81  balloon  launches  were 
made  in  September  of  1983,  43  from  Pt.  Mugu,  14  from  El  Segundo  and  24  from 
San  Nicholas  island.  The  rawinsonde  measurements  were  made  at  rou^ily  0.1 
km  increments.  The  pressure,  temperature  and  wind  velocities  were  then 
smoothed  using  a  cubic  spline  fit.  5  The  fitted  data  was  then  used  to 
evaluate  Cfi2  from  the  models  discussed  in  the  next  section. 

3.0  ANALYSIS 

The  model  used  for  the  free  atmosphere  was  that  developed  by  VanZandt  and 
coworkers.7  An  extensive  discussion  of  this  model  and  how  it  can  be  used 
in  conjunction  with  rawinsonde  data  is  provided  in  reference  6.  The  upper 
limit  used  was  arbitrarily  chosen  as  4  km  but  oould  easily  be  extended  up 
to  where  the  balloon  bursts,  which  is  typically  25  to  30  kilometers.  It  is 
assumed  that  the  lower  limit  on  this  model  is  at  the  bottom  of  the  first 
strong  temperature  inversion  layer.  This  should  mark  the  transition  from 
the  boundary  layer  to  the  free  atmosphere. 

The  model  used  for  the  boundary  layer  contribution  is  due  to  Kaimal  and 
coworkers.3  As  modified  fcy  Walters  and  Kunkel,®  the  model  takes  the  form: 

(Z/Z0)~V3  Z0f  ^  0>5Zi 

|  <-5Zi/Z0"V3  0*5zi  4  Z  <0.7Zi 

\ 

2.9(.5Zi/Z0)-V3(Z/Z.)3  0.7Zi  ^  Zi 


C„2(Z) 


!  9 


The  model  is  expected  to  work  best  when  there  is  a  well-defined  temperature 
inversion  present,  which  was  true  for  about  2/3  of  the  launches.  There  are 
two  parameters  built  into  this  model:  Zif  the  height  of  the  inversion 
layer,  and  Cn2(Z0),  the  value  of  Cn2  at  some  height  Z0. 

A  computer  algorithm  is  used  to  determine  froni  the  imputed  smoothed 
temperature  values  (at  0.1  km  increments)  by  searching  for  the  first 
occurrence  (if  any)  of  a  negative  temperature  gradient.  This  height  is 
taken  as  the  height  of  the  inversion  base  if: 

a)  the  temperature  gradient  remains  negative  for  the  next  two  0.1  km 
increments  and  drops  by  at  least  1°, 

or 

b)  the  temperature  drops  by  more  than  2°. 

When  a  temperature  inversion  is  found,  has  been  of  order  0.5  km,  with 
some  occurrences  of  up  to  order  2.5  km.  No  inversion  is  found  in  about 
1/3  of  the  launches. 


A  series  of  temperature  plots  were  done  for  three  launches  at  Pt.  ttigu  on 
September  28,  1983.  These  are  representative  of  a  well-defined,  relatively 
high-altitude,  inversion  layer.  This  layer  appears  to  move  upward  and  lose 
definition  as  the  day  progresses.  Another  plot  illustrates  a  launch  for 
which  there  is  apparently  no  inversion,  while  a  third  plot  illustrates  a 
weakly  defined  inversion  layer  at  about  0.7  km,  which  does  not  meet  the 
selecting  criteria  of  the  above  algorithm.  This  is  probably  all  right, 
since  this  inversion  is  so  weakly  defined  that  the  Kaimal  model  may  not  be 
applicable. 

The  second  parameter,  Cn2(Z0)  (i.e.,  an  accurate  measure  of  Cn2  dcwn  in  the 
stable  region  of  the  boundary  layer) ,  was  not  available.  In  the  absence  of 
such  a  measurement,  the  VanZandt  model  was  used  to  calculate  Cn2  at  the 
bottom  of  the  inversion  layer.  We  then  require  that  the  Kaimal  model  give 
the  same  value  at  this  point.  With  this  modification  we  have: 


cn2(z> 

cn2(z0 


.136  (Z/Zi)"V3 
.343 

wo 3 


Z  1  .5Zi 


20 


•5zi  i  .7Zi 
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where  Cn  (Z^  is  obtained  from  the  VanZandt  model.  Should  accurate  values 
of  Cfi2  at  an  appropriate  height  become  available,  we  could  easily  go  back 
to  the  original  formulation. 


4.0  RESULTS 

For  25  of  the  80  launches,  a  well-defined  temperature  inversion  was  located 
above  1  km.  For  each  of  these  we  calculated  Cn2  from  4  km  down  to  the 
inversion  layer  using  the  VanZandt  model.  We  then  forced  the  Kaimal  model 
to  match  the  VanZandt  model  at  this  point  and  then  used  it  to  calculate  Cn2 
down  to  0.1  km. 

Hie  upper  limit  of  4  km  is  quite  arbitrary  in  that,  typically,  the  data 
needed  for  the  VanZandt  model  is  available  up  to  at  least  20  km.  Hie  lower 
limit  of  0.1  km  is  somewhat  arbitrary  in  that  it  is  not  dear  how  far  down 
the  Kaimal  model  can  be  used.  The  use  of  relatively  high  (above  1  km) 
inversion  layers  was  made  to  show  the  contrast  between  the  boundary  layer 
and  the  free  atmosphere. 

For  contrast,  we  have  also  the  Kaimal  model  to  calculate  Cn2  from  the 
inversion  layer  down  to  0.1  km.  However,  it  must  be  remembered  that  this 
model  is  not  valid  in  this  regime. 

A  typical  result  for  Cn2  along  with  the  associated  temperature  profile 
shews  a  well-defined  inversion  beginning  at  2.5  km.  Graphically,  below  the 
inverison  layer,  the  improper  extension  of  the  VanZandt  model  is  noted  for 
comparison  to  the  better  suited  Kaimal  model. 

Analysis  of  three  nearly  simultaneous  launches  from  three  launch  sites  show 
that  for  each  set,  although  the  structure  of  the  temperature  profile 
through  the  inversion  differ,  the  height  is  approximately  the  same.  A  peak 
in  Cn2  is  also  common  for  each  of  the  three  launches. 

It  might  be  more  desirable  to  use  the  original  formulation  of  the  Kaimal 
model  instead  of  the  one  used,  even  though  the  two  are  mathematically 
equivalent.  Forcing  the  Kaimal  model  to  join  the  VanZandt  model  will 
affect  all  the  other  values  obtained  belcw  Z ^  p^r  a  iaunch  analyzed. 


changing  2^  to  2.9  km  or  3.1  km  will  affect  the  Cn2  values  below  Z^  by 
several  orders  of  magnitude,  were  an  accurate  value  of  Cn2  available  for  z 

<  this  problem  could  be  avoided. 

5.0  SUMMARY 

The  high  resolution  rawinsonde  data  obtained  at  the  three  stations  in  the 
Los  Angeles  area  provided  extremely  good  examples  of  well  established 
temperature  inversion  layers  which  generally  were  present  in  the  early 
afternoon.  This  fits  well  with  conditions  described  by  Kaimal 3  and  the 
application  of  the  modified  Kaimal  model  by  Walters  and  Kunkel8. 

We  have  calculated  CR2  as  a  function  of  altitude  when  a  well-defined 
temperature  inversion  layer  is  present  at  altitudes  of  1-3  km.  We  used  the 
VanZandt  model  down  to  the  inversion  layer  and  the  Kaimal  model  below.  One 
of  the  parameters  in  the  Kaimal  model  is  determined  from  the  VanZandt 
model;  the  other  parameter  (the  inversion  layer  height)  is  obtained  from 
the  smoothed  temperature  profile.  Should  independent  measurements  of  the 
first  of  these  parameters  become  available,  our  procedure  can  easily  be 
modified  to  take  advantage  of  it. 
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RAWINSCMDE  DATA  PROCESSING 

1.0  INTRODUCTION 

During  the  past  year,  BEDFORD  RESEARCH  personnel  transferred  raw insonde 
balloon  data,  called  "Cl earl"  and  "Clear2",  from  floppy  disks  to  the  CDC 
Cyber  System.  These  files  were  edited  and  then  damped  to  an  ANSI  standard 
labeled  tape  which  was  later  brought  to  White  Sands,  New  Mexico.  The 
following  is  an  in-depth  report  on  the  details  of  this  project. 


2.0  DATA  PROCESSING 

2.1  Qyeoiew 

The  original  eighty-eight  data  files  were  converted  to  ASCII  code  and  put 
on  five  and  one-quarter  inch  floppy  disks.  These  files  were  transferred  to 
the  Cyber  System  (see  Appendix  A  for  a  list  of  the  files) .  They  were  kept 
as  direct-access  files  because  the  total  prus  for  indirect-access  files  for 
that  account  was  insufficient.  The  following  describes  the  original  data 
in  detail. 

2.2  Pata_Descxiptiton 

Each  Cl earl  and  Clear2  data  file  consists  of  three  parts.  The  first  two 
are  header  records;  the  third  consists  of  the  actual  data  records.  These 
parts  are  described  below  in  more  depth. 

The  first  header  record  consists  of  five  fields.  These  fields  contain  the 
date,  time,  serial  number,  latitude,  and  longitude.  The  following  is  an 
example  of  the  first  header  records  (see  Appendix  B  for  the  format) : 

09-07 8421:11: 3 4nl 882  +32.40000+106.3700 


The  second  header  record  is  divided  into  two  records  for  writing  purposes. 
The  first  contains  these  surface  parameters:  altitude,  pressure,  temper¬ 
ature,  humidity,  wind  speed,  wind  direction,  C  n2,  and  10  o  ex2.  The  second 
contains  R  ^  ,  0^,  vrms,  zfinal,  grot  offset,  and  the  number  of  data 
records  (see  Appendix  C  for  formats) .  The  following  is  an  example  of  the 
second  header  record: 


1216.2  875.1+31.6  +  22.  6.2+300.+4.43E-14t0.00E+00 

+2.8  +3.4  +16.724.75  +6  999 

The  third  part  consists  of  the  actual  data  records.  In  the  original  Clear 
data,  the  fields  were  altitude,  pressure,  temperature,  humidity,  C^2,  the 
logarithm  of  the  average  Cn2,  a  of  the  logarithm  of  the  average  Cn2,  the 
east  displacement,  the  north  displacement,  wind  speed,  and  wind  direction 
(see  Appendix  D  for  the  formats  of  the  original  files).  The  following  is 
an  example  of  a  data  record  in  the  Clear 1  and  Clear 2  files: 

11.000  238.3-40.0  +21+2 . 27E-1 8+1 . 7 4E-1 8+1 . 27E+00  -9.99  -9.99  +23.1+271 


The  original  files  were  edited  by  a  FORTRAN  program.  The  new  files  have 
the  same  fields,  excluding  the  logarithm  of  the  average  Cn2,  a  of  the 
logarithm  of  the  average  Cn2,  the  east  displacement,  and  the  north 
displacement  (see  Appendix  E  for  a  list  of  the  new  files  and  Appendix  F  for 
the  formats  of  the  new  files) .  The  following  is  an  example  of  a  data 
record  in  the  Cl earl  and  Clear2  files: 

11.000  250.2-40.0  +21+2.27E-18  +23.1+271 

A  few  items  concerning  the  new  files  should  be  noted.  If  there  was  no  data 
gathered  for  the  humidity,  Cn2,  or  the  wind  direction,  the  field  received 
the  value  -999.  If  the  wind  speed  had  no  data  for  that  altitude,  the  value 
-999.9  was  given  to  the  data.  Also,  the  file  Om5315b  contains  two 
different  data  sets,  which  run  together  due  to  no  end-of-file  marker.  The 
original  data  was  in  this  format  for  this  file. 


2r> 


2.4  Data  Tape  Creation 


The  otigiral  and  new  files  were  dumped  to  CC  tapes.  Another  tape,  which 
was  taken  to  White  Sands,  New  Mexico,  has  the  new  files  on  it.  The  tape 
sent  to  White  Sands  was  an  ANSI  standard  labeled  tape.  It  was  a  nine-track 
tape  with  a  density  of  6250  bits  per  inch.  There  was  a  volume  header 
label,  a  header  label,  and  and  end-of-file  label  on  the  tape.  No 
end-of-volume  label  was  necessary. 

Each  file  was  copied  to  the  tape  in  coded  form  until  the  end-of-file  marker 
was  reached.  Each  end-of-file  marker  determines  the  end  of  one  file  and 
the  beginning  of  another.  However,  each  file  is  considered  one  record. 
The  maximum  block  size  was  set  to  10,000. 


3.0  SUMMARY 

The  transferral  of  the  Clearl  and  Clear2  files  to  the  Cyber  system  took 
place  in  May,  1985.  The  editing  and  creation  of  the  tape,  along  with  the 
documentation,  was  completed  in  June. 


The  tape  which  was  sent  to  White  Sands  may  be  recreated.  It  may  be  easier 
to  use  the  present  CC  tapes,  or  to  create  a  CC  tape  again. 
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UsJL_of_Files 


M0503 

M1888 

M0506 

Ml  889 

H0507 

M1890 

M0508 

M1891 

M0516 

Ml  892 

M0517 

M1893 

M0518 

Ml  895 

M0522 

Ml  896 

M0524 

M2694 

M0527 

M2703 

M0534 

M2710 

M0543 

M2711 

M0546 

M2712 

M0547 

M2754 

M1692 

M2757 

Ml  878 

M2758 

Ml  87  9 

M2759 

Ml  880 

M2762 

M1882 

M2763 

Ml  884 

M2764 

Ml  885 

M2766 

Ml  886 

M2769 

M2770 

M64L5 

M2771 

M64L6 

M6371 

M64L7 

M6372 

M6420 

M6373 

M6425 

M6381 

M6548 

M6383 

M6550 

M6384 

M6555 

M6385 

M6557 

M6386 

M6558 

M6387 

M6561 

M6388 

M6562 

M6389 

M6563 

M6392 

M6564 

M6395 

OM5308 

M6398 

OM5315b 

M6403 

OM5316 

M6406 

OM5317 

M6408 

OM5319 

M6409 

OM5321 

M64L0 

OM5322 

M64L2 

OM5323 
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APPENDIX  B 


First  Header  Record 


field 

columns 

format 

date 

1-8 

12.2, '-M2, '-M2 

+++ 

time 

10-16 

12.2, ,12.2, 1 ,12.2 

+++ 

serial  number 

17-21 

A5 

+ 

latitude 

26-33 

F8.4 

++ 

longitude 

34-42 

F9.4 

++ 

*  OH  files  have  six-character  serial  numbers,  and  OM5315A  has  a  seven- 
character  serial  number. 


+  Character  format  with  specified  length 
++  Single  precision  floating  point  real  without  exponent 
+++  integer  with  specified  length 


APPENDIX  C 


Second  Header  Record 


PartjL 

field  columns  format 


surface  altitude 

1-6 

F6.1 

** 

surface  pressure 

7-12 

F6.1 

** 

surface  temperature 

13-17 

F5.1 

** 

surface  humidity 

18-22 

F5.0 

** 

surface  windspeed 

23-26 

F4.1 

** 

surface  wind  direction 

27-31 

F5.0 

** 

Cn 2 

28-36 

E9.2 

* 

10  n  ex2 

37-45 

E9.2 

* 

Part_2 


R  ^ 

1-6 

F6.1 

** 

7-12 

F6.1 

** 

VRMS 

13-18 

F6.1 

** 

z final 

19-23 

F5.2 

** 

gmt  offset 

2  * ** ***-26 

13 

*** 

#  of  data  records 

27-30 

14 

*** 

*  Single  precisian  floating  point  real  with  exponent 

**  Single  precision  floating  point  real  without  exponent 

***  Integer  with  specified  length 
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Original  Data  Format 


field 

colunns 

format 

altitude 

1-6 

F6.3 

** 

pressure 

7-12 

F6.1 

** 

temperature 

13-17 

F5.1 

** 

humidity 

18-21 

14 

*** 

Cn2 

logarithm  of  average  Cn2 

22-30 

E9.2 

* 

31-39 

E9.2 

* 

a  of  logarithm  of  average  Cn2 

40-48 

E9-2 

* 

east  displacement 

49-55 

F7.2 

** 

north  displacement 

56-62 

F7.2 

** 

wind  speed 

63-68 

F6.1 

** 

wind  direction 

69-72 

14 

*** 

*  Single  precision  floating  point  real  with  exponent 
**  Single  precision  floating  point  real  without  exponent 
***  Integer  with  specified  length 
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M0503B 

M0506B 

M0507B 

M0508B 

M0516B 

M0517B 

M0518B 

M0522B 

M05243 

M0527B 

M053® 

M0543B 

M0546B 

M0  547B 

M1692B 

M1878B 

M1879B 

M1880B 

M1882B 

MI884B 

M1885B 

M1886B 


r.iat  of  Files 


M1888B 

M1689B 

M1890B 

M189IB 

M1892B 

M1893B 

M1895B 

MI896B 

M2694B 

M2703B 

M2710B 

M2711B 

M2712B 

M2754B 

M2757B 

M2758B 

M2759B 

M2762B 

M2763B 

M2764B 

M2766B 

M2769B 


M2770B 

M277IB 

M6371B 

M6372B 

M6373B 

M6381B 

M6383B 

M6384B 

M6385B 

M6386B 

M6387B 

H6388B 

M6389B 

M6392B 

M6395B 

M6398B 

M6403B 

M6406B 

H6408B 

M6409B 

M64L0B 

M64L2B 


M64L5B 

M6416B 

M6417B 

M6420B 

M6425B 

M6548B 

M6550B 

M6555B 

M6tiS7B 

M6558B 

M6561B 

M6562B 

M6563B 

M6564B 

^M5308B 

OM5315B 

OM5316B 

OM5317B 

OM5319B 

OM5321B 

OM5322B 

OM5323B 
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New  Data  Format 


field 

columns 

format 

altitude 

1-6 

F6.3  ** 

pressure 

7-12 

F6.1  ** 

temperature 

13-17 

F5.1  ** 

humidity 

18-21 

I4  *** 

C  2 

22-30 

E9.2  * 

n 

wind  speed 

31-36 

F6.1  *** 

md  direction 

37-40 

14  *** 

*  Single  precision  floating  point  real  with  exponent 

**  Single  precision  floating  point  real  without  exponent 

Integer  with  specified  length 


*** 


GAUSSIAN  SMOOTHING  OF  C2  OBSERVATIONS 
FROM  DAY  AND  NIGHT  FLIGHTS 


1.0  INTRODUCTION 

1.1  Q 

This  report  outlines  the  ADP  support  that  Bedford  Research  Associates 
provided  the  Atmospheric  Optics  Branch  (OPA)  of  the  Air  Force  Geophysics 
Laboratory  (AF3L) .  It  describes  the  nature  of  the  analyses  as  well  as  the 
software  used  in  this  investigation. 

1.2  Background  and  Objectives 


One  measure  of  atmospheric  turbulence  is  the  optical  refractory  constant, 
Cn2.  This  index  is  subject  to  a  variety  of  meteorological  events  such  as, 
temperature,  wind,  pressure,  and  relative  humidity.  Each  of  these  events 
have  been  used  in  the  prediction  of  Cn2  variation.  One  condition  that  has 
received  little  empirical  study  is  the  differential  effects  of  day  and 
night  sampling  on  Cn2  measures.  It  is  the  purpose  of  this  investigation  to 
examine  and  report  the  daily  and  nightly  patterns  of  Cp2  data. 

The  tasks  developed  to  analyze  these  data  are  presented  below: 


*  Download  data  files  from  tape  to  disk 

*  Edit  downloaded  files  to  exclude  cases  with  missing  data  and  adjust 
noise  C2  values. 

*  Write  the  algorithms  to  calculate: 

-  Running  means  of  Cn2  for  individual  flights 

-  Grand  (pooled)  C  2  means  from  the  individual  flights 

-  Gaussian  smoothed  Cn2  scores  from  individual  and  pooled  data  sets 

-  Regression  coefficients  for  the  pooled  mean  values  using  log 
quadratic  and  nonlinear  equations 

*  Plot  individual  and  pooled  CQ2  values  by  altitude 

*  Write  an  output  file  containing  altitude  and  gaussian  smoothed 
scores  to  tape 

1.3  fiepoct  Outline 


This  report  is  divided  into  three  sections.  The  first  section  provides  a 
description  of  the  software  designed  to  perform  the  tasks  listed  in  Section 
1.2.  In  this  part  of  the  report,  the  software  requirements,  functions 


executed  by  the  main  program,  and  structure  of  the  original  and  edited  data 
files  are  reviewed.  The  next  section  identifies  and  describes  the  results 
of  the  analyses  performed.  The  third,  and  final  section,  summarizes  the 
findings  of  this  study. 

2.0  SOFTWARE  SYSTEM 

2.1  flyerviw 

This  section  outlines  the  functional  characteristics  of  the  software 
developed  to  compute  the  running  means  and  gaussian  smoothed  Cn2  values  from 
individual  and  aggregated  day  or  night  flights.  It  describes  the  software 
requirements  used  to  create  this  system.  It  also  delineates  the  operations 
(i.e.,  inputs,  processing,  and  outputs)  within  the  main  program  and  presents 
the  format  of  the  data  files. 

2.2  Requirements 

In  the  initial  stages  of  software  development,  a  number  of  programming 
requirements  were  specified.  These  requirements  were  used  in  constructing 
the  algorithms  and  source  code  for  the  main  program.  Each  of  these 
constraints  are  discussed  below: 

*  Records  containing  missing  altitude  and/or  Cn2  data  were  omitted  from 
any  computations.  These  cases  provided  inaccurate  information  that 
could  distort  the  calculations  of  mean  or  smoothed  Cn2  scores. 

*  Cn2  observations  with  noise  were  adjusted  using  the  following 
equation: 

Cn2  =  (79.9  *  P/0T2)  *  10-*) 2  *  Oj.2 

where  P  =  pressure  (mb) 

T  =  temperature  (K) 

Cj2  =  (.002)2 

This  formula  resulted  in  values  of  Cn2  that  could  be  used  in  subsequent 
computations.  If  the  pressure  and/or  temperature  were  missing,  the  case  was 
excluded  from  further  processing. 

*  Mean  and  gaussian  values  were  estimated  on  cases  that  contained  an 
altitude  window  (or  standard  deviation)  of  20  or  40  meters  (m) . 
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Thus,  to  be  included  in  these  calculations,  an  altitude  level  of  10 
or  20  m  above  the  lowest,  and  below  the  highest,  points  was  required. 
Cases  that  did  not  fall  within  this  region  were  omitted  because  they 
did  not  contain  a  full  window. 

*  The  gaussian  jpoothing  of  individual  day  and  night  flights  were  done 
on  the  raw  Cn2  scores  contained  within  ea'n  file.  However,  the  grand 
smoothing  of  flights  was  performed  on  the  average  Cn2  values. 

*  Results  were  plotted  on  the  Calcomp  Plotter  or  written  to  a  tape  file 
containing  the  pooled  smoothed  cJ2  scores  plus  their  corresponding 
altitude  levels  for  day  and  night  flights. 

2.3  Software  Functions 

All  software  was  developed  and  executed  on  the  APGL  Control  Data  Corporation 
(CDC)  170/750  Cyber  mainframe  running  under  NOS  2.3.  The  source  code  for 
the  main  program  was  written  in  ANSI  standard  Fortran  (Version  5.0)  . 
Routines  were  also  coded  using  NOS  command  language  for  file  management  and 
IMSL  libraries  were  used  to  obtain  the  regression  coefficients. 

The  software  operations  were  implemented  on  the  individual  and  pooled  data 
files.  Therefore,  the  following  discussion  will  address  each  of  these 
activities  independently.  A  flowchart  is  presented  in  Figure  2.3.1  to 
describe  the  ADP  activity. 

2.3.1  individual  Flight  Operations 

2. 3. 1.1  Inputs 

The  inputs  to  the  individual  flight  computations  were  obtained  from  the 
downloaded  tape  files.  These  files  were  in  binary  format  and  are  described 
in  more  detail  in  Section  2.4.  The  inputs  to  the  algorithm  were  the  flight 
number,  altitude,  temperature,  pressure,  and  raw  Cn2  observations. 

2. 3.1.2  Processing 

The  data  was  first  sorted  by  altitude  in  ascending  order  and  then  read  into 

an  algor ithm  that  calculated  mean  and  smoothed  Cn2  scores.  At  the  beginning 

of  this  process,  the  integrity  of  each  observation  was  checked.  Records 

with  missing  or  outlying  values  were  deleted  from  the  file  and  those  with 

noise  Cn2  recording  were  corrected.  The  records  that  passed  through  this 

editing  stage  were  entered  into  an  array  of  valid  cases.  An  altitude  window 

was  constructed  for  each  element  in  the  array.  If  data  points  contained  a 

full  window  of  20  or  40  m,  they  were  used  to  compute  the  mean  and  smoothed 

C  2  values, 
n 
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2. 3.1. 3  Output 


The  results  of  individual  flight  operations  were  either  plotted  or  stored  on 
disk  for  further  analysis.  Plots  were  generated  for  the  3  day  and  6  night 
flights  first  using  a  window  of  20  m  and  then  for  40  m.  Thus,  a  total  of  18 
plots  were  created  from  the  individual  flight  data. 

Files  containing  the  flight  numbers,  altitudes,  and  mean  and  smoothed  Cn2 
scores  were  also  produced  at  this  stage  of  the  analysis.  These  files  were 
then  concatenated  into  one  day  and  one  night  file  fran  which  pooled  mean  and 
smoothed  Cn2  measures  were  constructed. 

2.3.2  Pooled  Flight  Operations 

2. 3. 2.1  inputs 

The  concatenated  files  of  day  and  night  flights  were  used  as  inputs  to  the 
pooled  flight  program.  The  file  description  is  presented  in  Section  2.4. 

2. 3.2.2  Processing 

The  concatenated  files  were  sorted  by  altitude  in  ascending  order  and  then 
passed  through  the  same  data  checks  that  were  performed  on  the  individual 
data  files.  Once  the  integrity  of  the  observations  he?  been  tested,  they 
were  passed  into  an  algorithm  that  computed  the  pooled  mean  Cn2  value  from 
the  individual  mean  scores.  The  pooled  mean  values  were  stored  and  used  as 
inputs  for  the  next  algorithm  that  computed  a  gaussian  smoothing  of  these 
values.  The  pooled  mean  Cn2  values  were  also  entered  into  a  log  quadratic 
and  nonlinear  equation. 

The  mean  and  smoothed  Cr2  values  were  created  twice.  The  first  time  a 
standard  deviation  of  20  was  used  and  the  second  time  a  standard  deviation 
of  40  was  employed. 

2. 3.2. 3  Output 

The  products  of  the  pooled  processing  were  either  plotted  or  stored  on 
disk/tape  files.  Three  plots  were  generated  with  altitude  assigned  to  the 
Y-axis  and  the  corresponding  Cn2  value  on  the  X-axis.  The  Cn2  treasure  was 
derived  from  either  the  gaussian  smoothing,  log  quadratic,  or  non-linear 
computations.  These  plots  were  each  created  from  the  files  computed  with  a 
window  of  20  or  40. 


37 


(Xitput  files  are  written  at  the  end  of  data  processing.  The  two  files 
generated  contained  the  altitude  and  pooled  mean  and  gaussian  smoothed  Cr2 
values.  The  only  difference  between  these  two  files  is  that  the  scores  were 
created  using  differ  ait  standard  deviations. 

2.4  Data  File  Characteristics 

2.4.1  Magnetic  Tape  Files 

Three  day  and  six  night  files  were  used  in  this  study.  These  files  resided 
on  an  ANSI-labeled  tape  at  1600  bpi.  The  tape  contained  22  files  written  in 
binary  format.  The  structure  of  these  files  is  presented  in  Appendix  A. 
All  files  had  the  same  record  format.  Each  file  consisted  of  five  header 
records  that  indentified  the  flight  parameters  (e.g. ,  launch  date  and  time) 
and  data  records  that  had  the  observations  for  each  altitude  level  sampled. 

2.4.2  Disk  Files 

Since  not  all  of  the  data  on  the  tape  files  were  necessary  for  the  present 
purposes,  records  were  truncated  when  written  to  disk.  Disk  files  cmly 
contained  the  flight  number,  altitude,  temperature,  pressure,  and  Cn2 
scores.  However,  the  pooled  data  sets  only  included  the  altitude,  mean  and 
gaussian  smoothed  Cn2  scores. 

3.0  DATA  ANALYSIS 
3*1  Overview 

This  section  reviews  the  results  of  the  analyses  conducted  within  this 
investigation.  It  first  describes  the  mathematical  approach  used  and  then 
reviews  the  tasks  completed. 

Quantitative  Methods 

Ttoere  were  two  procedures  used  to  summarize  the  observations  in  the  data 
files.  The  first  one  was  to  obtain  a  mean  Cn2  scores  for  a  given  window  as 
defined  below: 
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where  Zj  represents  the  altitude  of  Cn2(zj),  Zi  is  the  altitude  of  the  mean 

cn2(zi)  score,  and  n  equals  the  number  of  scores  that  fall  within  the 
window.  This  sample  mean  denotes  an  unbiased  estimate  of  the  population 
mean. 

The  second  procedure  used  to  aggregate  the  data  was  a  gaussian  smoothing 
technique.  This  approach  utilized  a  least-squares  estimate  and  is  expressed 
as  follows: 

A(Zi>zj)Cn2(zi) 

A<Zi,z;j) 

where: 

A(ZirZj)  =  1/-  2  e-l/2(Z*‘~zi) 2/ 

The  mean  Cn2  values  derived  from  the  first  equation  were  used  as  inputs  to 
the  log  quadratic  and  non-linear  equation.  The  log  quadratic  formula  is 
shown  below: 

109  Cn2(zi)  =  A  +  BZi  +  CZi2 

where  A,  B,  and  C  are  coefficients  of  the  quadratic  and  Z^  is  the  altitude 
of  case  i.  Zi  ^  2 j2  are  assumed  to  be  orthogonal. 

The  nonlinear  ( Newton- Raphson)  equation  used  in  this  investigation  is 
represented  as: 


Log  Cn2(Zi)  a  a  +  BZi  +  CZ^2  +  De*P  (-l/2(Zi-fy'fr) ) 2 


where  A,  B,  C,  D,  E,  and  F  are  coefficients  within  the  expression. 


3.3  Results 

3.3.1  Individual  Flights 

The  product  of  the  gauss ian  smoothing  of  an  individual  day  flight  are 
presented  in  Figures  3. 3. 1.1  and  3. 3. 1.2.  Both  figures  display  the 
logarithm  of  smoothed  Cn2  scores  on  the  X-axis  and  their  corresponding 
altitude  levels  (in  tan)  on  the  Y-axis.  The  Cn2  values  for  the  first  figure 
were  computed  using  an  altitude  standard  deviation  of  20  m  while  the  second 
figure  used  a  standard  deviation  of  40  m.  A  series  of  plots  were  produced 
for  different  flights.  Flight  2X0006  demonstrated  the  most  Cn2  variation. 
Incrementing  the  window  from  20  to  40  m  shows  a  pronounced  effect  on  the 
profile  of  smoothed  C n2  values. 

The  gaussian  smoothing  of  an  individual  night  file  are  shewn  in  Figures 

3. 3. 1.3  and  3. 3. 1.4.  Once  again  the  smoothed  Cn2  scores  were  plotted  on  the 
X-axis  and  the  altitude  on  the  Y-axis  and  the  smoothed  Cn2  values  were 
computed  using  either  a  20  or  40  ro  window.  The  differences  between  the  20 
and  40  m  computations  are  noticeable.  The  40  m  window  displays  less 
variation  in  Cn2  than  the  20  m  window. 

3.3.2  Fooled  Flights 

The  outputs  resulting  from  the  operations  performed  on  the  concatenated  day 
flights  are  presented  in  Figures  3.3. 2.1.  Figure  3. 3. 2.1  displays  a  graph 
of  the  pooled  mean  Cn2  scores  by  their  corresponding  altitudes.  The  graph 
was  computed  using  the  standard  deviation  of  20  m.  Figure  3. 3.2.2  shews  the 
results  of  the  gaussian  smoothing  of  pooled  mean  Cn2  night  flights. 

Figure  3. 3. 2. 2  depict  the  pooled  mean  and  gaussian  smoothing  of  the  pooled 
mean  Cn2  for  the  concatenated  night  flights  having  a  20  m  window.  From 
inspection  of  the  20  m  and  40  m  plots,  the  decision  was  made  to  only  include 
the  20  m  windowed  data  set  in  the  quadratic  and  nonlinear  fits. 

The  coefficient  used  to  plot  the  quadratic  and  nonlinear  fits  of  pooled  mean 
Cn2  values  are  shown  in  Table  3. 3.2.1.  Using  the  coefficients  displayed  in 
the  top  portion  of  this  table,  the  fitted  lines  were  estinated  for  the  day 
and  night  flights  and  plotted  in  Figure  3. 3. 2. 3  along  with  the  smoothed  Cn2 
lines.  Similar  plots  were  drawn  using  the  coefficients  shewn  at  the  bottom 
part  of  Table  3. 3. 2.1  were  drawn  in  Figures  3. 3.2. 4  and  3. 3. 2. 5.  These 
plots  represent  the  nonlinear  fit  of  day  and  night  data,  respectively. 
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Table  3. 3.2.1 

Table  of  Coefficients  Used  in 
Quadratic  and  Non-Linear  Fits 


Coefficients 


Eit 

Label 

u ax 

Night 

Quadratic: 

A 

.2709820416  40  81E-15 

.66 40331056 8815E-16 

B 

-.2934551 9050923E-16 

-. 13270127 983623E-17 

C 

.789367  447831 21E-18 

.77160351 583 446E-19 

Nonlinear: 

A 

-.773795E+2 

-. 76380 4E+2 

B 

.152222 

-.434465E-2 

C 

.128323E-2 

-.418050E-2 

D 

.544582E+5 

.528689E+5 

E 

•199407E+5 

.195242E+5 

F 

•  54L586E+4 

.5301 45E+4 
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APPENDIX  A 


Data  Format  of  Mag  Tape  File: 
VSN  =  CC4459 


Word  Nypiter 


Content 

iieadet  jhq  l: 

Name 

Launch  date 
Launch  time 
Reference 
Minimum  Frequency 
Maximum  Frequency 
Run  Date 


Format 


A10 

A10 

F10.2 

F10.2 

F10.2 

F10.2 

A10 


R0 

*30 

T33  4. 

Tgwet 

^ry 

Pa 


h 

Run  Date 
*  Words/Record 


F10.2 

F10.2 

F10.2 

F10.2 

F10.2 

F10.2 

F10.2 

F10.2 

A10 

110 


Thermosonde 


#  words  following  low  gain  14 
1st  time  division  F10.2 

1st  thermosonde  correlation  F10.2 

2nd  time  division  F10.2 

2nd  thermosonde  correlation  F13.2 


—  Repeats  for  #  words  specified  in  word  1 


#  words  following  high  gain  14 
1st  time  division  14 

1st  thermosonde  correlation  F10.2 

2nd  time  division  F10.2 

2nd  thermosonde  correlation  F10.2 


Vi*] 


Record  Number 
Time 

Frequency  Reference 
CINT 

Temperature  (C°) 

C.  Temperature 
CH 
C  CH 

.Contents: 

Altitude 

C.  Altitude 

Frequency  pressure 

KLASS 

IPRESS 

Pressure 

Frequency  -  (Thermo-low) 
Thjjrmosonde  CT 

Frequency  (Thermo-high) 
Thermo  CT 


S’2  combined 
[titude 
Position  East 
Position  North 
Vel.  East 
Vel.  North 
Vel.  Vertical 
Vel.  Horizontal 
SHOR 

Smoothed  Cvj2 
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TURBOPAUSE  HEIGHT  STUDY 


1.0  INTRODUCTION 

The  purpose  of  this  study  was  to  determine  if  a  relationship  exists 
between  turbopause  height  and  solar  flux  (measured  by  F^0  7) . 


2.0  DATA  PROCESSING 

2.1  Data  Bsssuftian 

Fifteen  years  of  rocket  data  were  collected  from  two  stations:  Eglin 
AFB  and  Heiss  Island.  Variables  recorded  included  date,  local  time,  and 

turbopause  height.  Apf  F10.7  measurements  have  been  integrated 

into  the  data  set  using  the  time  and  data  information. 

2.2  Data  Analysis 

2.2.1  Overview 

Raw  data  was  inspected  for  outliers  using  scatter  diagrams  of  turbopause 
height  and  other  variables.  Some  data  records  were  removed  if  inaccurate, 
or  alternatively  the  turbopause  height  was  adjusted.  Kp  values  were  also 
corrected  as  necessary.  A  final  collection  of  70  records  were 
deemed  valid  for  further  consideration. 

The  data  was  split  into  seasons,  defined  as  follows: 

Winter  —  Dec  1  -  Mar  15  26  records 

Spring/Sunmer  —  Iter  16  -  Sep  15  25  records 

Fall  —  Sep  16  -  Nov  30  19  records 

2.2.2  Fall  Data 

The  fall  data  was  removed  from  further  analysis  due  to  its  lack  of 
significance  in  determining  a  relationship  between  solar  flux  and 
turbopause  height.  This  was  determined  by  the  Initiator  from  previous  work 
on  the  fau.1  data. 


Initial  descriptive  statistics  (means,  variances,  and  correlation  matrices) 
of  the  winter  and  spring/summer  data  sets  were  dene  on  all  the  variables. 
For  the  sake  of  completeness,  the  fall  data  was  included  in  this  prelimin¬ 
ary  analysis. 

2.2.3  Winter  Data 

A  descriptive  statistics  analysis  was  done  to  examine  the  characteristics 
of  the  winter  data  set.  Significant  correlation  coefficients  were  observed 
in  the  correlation  matrix.  They  are  tabulated  as  follows: 

Turbopause  Height  with  Kp^  .4157 

Ap3  with  Kpl  17731 

Kp2  and  -.4569 

A  stepwise  multiple  regression  of  the  winter  data  was  done.  The  Kpi 
variable  entered  the  regression  equation  at  a  .05  level  of  significance. 
The  resulting  equation: 

Turb.  Ht.  =  105.42  +  1.5(Kpi) 

has  a  multiple  correlation  coefficient  of  .173,  which  is  not  a  highly 
significant  model  result. 

2.2.4  Spring/Suimver  Data 

The  original  data  set  of  spring/summer  was  examined.  The  descriptive 
statistics  along  with  a  set  of  scatter  diagrams  are  presented.  There  were 
no  significant  relationships  present  between  any  of  the  variables  in  the 
correlation  matrix  of  the  spring/summer  data  set.  A  stepwise  multiple 
regression  analysis  was  performed  and,  likewise,  no  variables  were  entered 
into  the  model  to  estimate  the  turbopause  height, 

2.2.5  Diurnally  Corrected  Data 

For  the  winter  and  spring/summer  data  sets,  a  diurnal  component  was 
detected  after  examination  of  local  time  vs.  turbopause  height  plots.  The 
diurnal  variation  of  each  season  was  removed  and  a  summary  statistics 
analysis  was  performed. 


Similar  results  to  those  in  2.2.3  were  found  in  the  diurnally  corrected 

winter  data.  That  is,  had  a  high  enough  correlation  (at  the  .05  level) 
to  be  included  in  a  stepwise  regression  model.  The  correlation 
coefficient  between  turbopause  height  and  was  .3889.  The  resulting 
equation: 

Turb.  Ht.  =  -2.19  +  1.22(Kpi) 

has  a  irultiple  correlation  coefficient  of  .151,  which  is  not  a  highly 
significant  model  result. 

The  diurnally  corrected  spring/summer  data  had  no  statistically  significant 
characteristics  for  accessing  the  relationship  between  solar  flux  and 
turbopause  height,  but  the  scatter  plot  of  F10>?  ^  turbopause  height 
revealed  an  interesting  structure. 

2.2.6  Yearly  Averaged  Data 

Yearly  turbopause  height  values  were  averaged  for  each  season.  Descriptive 
statistics  and  scatter  diagrams  were  utilized  to  determine  if  any  relation¬ 
ship  existed  between  turbopause  height  and  solar  flux  (F10  7) .  &  hand 

drawn  curve  was  placed  over  the  turbopause  height  vs  F10  7  graph  to  examine 
the  trend  that  turbopause  height  had  with  time.  A  baric  pattern  developed: 
in  winter  a  counter  clockwise  rotation  of  turbopause  height  was  seen 
through  time,  .and  Id  the  spring/ summer  o  general  clockwise  rotation 
was  seen  through  time* 

2.2.7  Station  Split  Data 

The  spring/summer  and  winter  data  vas  split  into  each  station:  Eglin  AFB 
and  Heiss  Island.  A  similar  hand  drawn  curve  was  overlaid  onto  the  turbo- 

pause  height  vs  F^y  graphs  to  examine  the  flew  of  turbopause  height  in 
relation  to  solar  flux  through  time.  Due  to  the  scarcity  of  data  in  each 
station's  season,  (n=9,12,13,or  17)  a  vague  ellipse  is  evident 
showing: 

winter  -  counter  clockwise  .rotation  of  tur^opause  height  in 
relation  to  solar  £iu*  through  time 

Spring/Summer  -  clockwise  rotation  of  turbopause  height  in 
relotion  to  solar  JLlua  through  rime* 


3.0  SUMMARY 

A  basic  relationship  between  turbopause  height  and  the  solar  rotational 
effect  has  been  determined.  Even  though  the  data  set  had  been  scarcely  and 
inconsistently  gathered  throughout  time  (70  valid  records  over  15  years) ,  a 
distinct  trend  in  the  direction  of  turbopause  height  in  relation  to  solar 
flux  through  time  was  evident.  It  is  understood  that  there  was  great 
expense  in  gathering  each  record,  which  consisted  of  a  rocket  launch. 


The  relationship  between  turbopause  height  and  the  solar  rotation  cannot  be 
described  by  a  formal  statistical  model  due  to  the  nature  of  the  data  set 
as  mentioned  in  the  previous  paragragh.  However,  again  stated,  the 
following  relationship  did  appear  founded.  Given  diurnally  corrected  and 
yearly  averaged  turbopause  height  values,  as  well  as  a  solar  flux  indicator 

(pi0.7) ,  the  following  relationships  were  noted: 


(1)  There  is  evidence  of  £  counter-clockwise  rotation  Of  the 
turbopause  height  in  relation  to  the  solar  rotation  during  the 
winter  (defined  from  December  1st  throucji  March  15th)  season. 


(2)  Similarly,  there  is  evidence  of  £  clockwise  rotation  Of  the 
turbopause  height  during  in  relation  to  the  solar  rotation  during 
the  spr ing/fiiut»ner  (defined  from  March  16th  to  September  15th) 
season. 


It  is  evident  that  more  data  must  be  gathered  in  order  to  substantiate  this 
claim. 
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INSTANTANEOUS  RAINFALL  ANALYSIS 

1.0  INTRODUCTION 
1.1  Overview 

This  report  summarizes  the  technical  and  analytical  support  provided  by 
Bedford  Research  Associates  (BRA)  to  the  Design  Climatology  Program  which 
is  directed  by  the  Tropospheric  Structure  Branch  of  the  Air  Force 
Geophysics  Laboratory  (APGL)  located  at  Han scorn  Air  Force  Base,  Bedford, 
MA.  It  describes  the  ADP  software  and  analyses  used  to  evaluate  the 
parameters  of  intense  rainfall. 


The  Design  Climatology  Program  develops  techniques  for  incorporating 
climatic  information  into  the  design,  operation,  and  evaluation  of 
aerospace  equipment.  One  method  under  current  investigation  is  the 
statistical  modeling  of  1-minute  rainfall  rates.  The  study  of  rainfall  is 
significant  because  rain  can  disrupt  normal  Air  Force  operations  in  many 
ways.  For  example,  heavy  intensities  of  precipitation  attenuate  the 
microwave  signals  transmitted  by  radar  systems  for  satellite  detection  and 
tracking,  surveillance,  and  reconnaissance.  Rainfall  also  adversely 
effects  the  operation  of  aircraft  through  erosion  of  leading  edges  (e.g. , 
helicoptor  blades)  or  malfunction  of  jet  engines.  Therefore,  a  better 
understanding  of  the  characteristics  of  rainfall  rates  will  enhance  the 
design  and  function  of  Air  Force  systems. 

Despite  the  significance  of  this  atmospheric  event  on  aerospace  systems, 
there  exists  little  empirical  data  on  the  characteristics  of  1-minute, 
instantaneous  rainfall  rates.  Tattelman  and  Grantham  (1982)  have  reviewed 
the  literature  and  concluded  that  estimates  of  instantaneous  rainfall  are 
restricted  by  inherent  methodological  weaknesses.  Specifically,  errors  in 
data  collection,  caused  by  inaccurate  measuring  devices  and  large  annual 
variations  in  the  rainfall  distribution  raise  questions  about  the  validity 
of  previous  findings.  The  only  technique  for  developing  an  accurate  model 
of  instantaneous  rainfall  is  to  compile  a  data  base  containing  continuous 
and  reliable  measurements  for  at  least  5  years.  This  strategy  should 
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compensate  for  the  experimental  error  resulting  from  poor  instrumentation 
and  yearly  variations. 

1.3  Objectives 

Hie  primary  objective  of  this  investigation  was  to  evaluate  the  frequency 
and  duration  of  1-minute  rainfall  rates  compiled  over  a  10  year  period. 
This  study  assessed  the  annual,  monthly,  and  seasonal  patterns  of  rain  in 
order  to  develop  a  more  comprehensive  understanding  of  its  properties.  More 
specifically,  the  tasks  conducted  under  this  project  were  to: 


*  Examine  the  contents  of  the  tape  files  to  determine  the  integrity  of 
the  data. 

*  Identify  rainfall  rates  ranging  between  .001  and  10.5  nn/min  and 
count  the  frequency  of  occurrence  at  various  intervals  and  durations. 

*  Aggregate  the  rainfall  by  annual,  monthly,  or  seasonal  groups  and 
plot  the  frequency  distribution  for  various  durations. 

*  Investigate  the  number  of  times  that  the  rainfall  duration  crossed 
seasonal  boundries. 

*  For  each  season  within  the  10  year  span,  determine  the  longest  dura¬ 
tion  of  rainfall  with  rates  between  .1  and  2.5  nnv/min  and  report  the 
year  and  length  of  that  duration  by  season  and  threshold  value. 

*  Repeat  the  above  analysis  for  each  of  the  months  within  the  data 
file  and  report  the  results  by  months  instead  of  seasons. 

*  Calculate  and  plot  the  Poisson  probability  distribution  of  seasonal 
rainfall  for  varying  duration  and  thresholds. 

*  Explore  the  relationship  between  the  rainfall  intensities  and  dura¬ 
tion  lengths  and  tabulate  the  findings. 

*  Tabulate  the  rainfall  rates  at  thresholds  ranging  from  .001  to  2.5 
nr/min  by  the  season  and  hour  of  the  day. 

*  Group  the  data  by  4-hour  intervals  and  plot  the  frequency  distri¬ 
bution  for  different  thresholds  and  seasons. 

*  Refine  the  output  generated  from  the  above  analyses  for  final 
delivery. 


The  ultimate  goals  were  to  establish  a  more  comprehensive  explanation  of 
the  properties  of  instantaneous  rainfall  and  to  develop  software  that  could 
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be  used  to  replicate  these  analyses  on  other  rainfall  data  bases. 

1.4  lfett»dology 

This  project  was  partitioned  into  four  distinct  phases.  The  first  phase 
involved  determining  the  structure  and  content  of  the  data  sets.  All 
relevant  magnetic  tapes  and  documentation  were  reviewed  to  insure  the 
integrity  of  the  elements  with  the  file.  Programs  were  also  run  to 
identify  the  file  characteristics.  Once  the  data  was  verified,  the  second 
phase,  establishing  software  requirements,  commenced.  During  this  process, 
the  specific  statistical  analyses  were  identified.  The  data  to  be  analyzed 
(e.g. ,  intensity,  duration)  and  the  form  of  the  output  (e.g.,  tables  or 
plots)  were  defined.  Based  on  these  specifications,  the  routines  were 
written  and  placed  within  the  main  program.  The  final  phase  involved 
reviewing  the  results  to  determine  the  subsequent  analyses. 

1.5  Esppcfc  .Organization 

The  remainder  of  this  report  describes  the  software  development  and  data 
analyses  used  in  the  rainfall  project.  In  the  next  section,  the  main 
program  and  data  files  are  presented.  The  operations  performed  by  the  main 
program  and  the  characteristics  of  the  rainfall  files  are  discussed.  The 
subsequent  section  explains  the  data  analyses  performed  under  this  project. 
The  final  section  provides  a  summary  of  the  results  and  considerations  for 
future  research  in  the  area  of  instantaneous  rainfall. 

2.0  SOFTWARE  SYSTEM 
2.1  Overview 

This  section  presents  a  functional  description  of  the  software  written  to 
evaluate  the  rainfall  data  files.  It  identifies  the  inputs,  processing, 
and  outputs  of  the  main  program  used  in  this  investigation.  In  addition, 
it  outlines  the  format  and  structure  of  the  data  files.  The  purpose  of 
this  section  is  to  identify  the  major  ADP  components  of  this  project. 
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2.2.1  Hardware 

All  ADP  work  was  performed  on  the  APGL  Control  Data  Corporation  (CDC)  Cyber 
mainframe.  This  system  had  a  170/750  processor  with  26 2K  words  of  primary 
memory.  It  accessed  14  disk  drives  for  a  total  storage  capacity  of 
approximately  ,5M  words.  Also,  there  were  11  tape  drives  (mixed  7-  and  9- 
track)  that  handled  tape  densities  between  556  and  6250  bpis.  Pour  line 
printers  (3  1200  1pm  or  1  Xerox  page  printer)  and  a  Calcomp  1055  drum 
plotter  were  available  to  generate  output  from  the  system. 

2.2.2  Software 

The  Cyber  mainframe  ran  under  the  NOS  2.3  operating  system.  It  supported  a 
variety  of  software  needed  for  text  editing,  file  manipulation,  and  data 
analysis.  The  Interactive  Facility  (IAF)  was  used  to  access  the  XEDIT  text 
editor  that  created  and  modified  program  files.  The  source  code  for  the 
main  program  was  written  in  ANSI  standard  Fortran  (Extended  Version  4.7  or 
5.0) .  Programs  were  also  written  in  NOS  command  language  to  perform  file 
maintainance.  Fortran  and  AFGL  Plotting  libraries  were  used  to  plot  the 
results. 


The  main  program  developed  for  the  rainfall  analysis  was  written  in  Fortran 
5.0  and  stored  on-line  in  a  indirect  access  file  and  off-line  on  a  backup 
tape.  The  remainder  of  this  section  explains  the  operations  performed  by 
the  main  program. 

2.3.1  Inputs 


Inputs  were  obtained  from  two  sources.  Hie  first  source  was  from  within 
the  batch  file.  These  inputs  defined  the  parameters  of  that  run.  These 
parameters  included  identification  of  the  plotting  profile  (linear,  semi- 
linear,  or  logarithmic),  definition  of  the  X-axis  origin  for  plotting, 
determination  of  the  number  of  duration  intervals  and  time  span  (e.g., 


seasonal,  monthly,  or  annual)  used  in  the  computation. 

The  second  source  of  input  was  the  1-minute  rainfall  data  files.  TVo 
direct  access  files  contained  rainfall  measures  sampled  at  1-minute 
intervals  for  a  combined  total  of  over  10  years.  These  two  files  were 
sequentially  attached  to  the  main  routine  each  time  it  was  run.  (A  more 
detailed  description  of  the  contents  of  these  files  is  presented  in  Section 
2.4) . 

2.3.2  Processing 

Meetings  were  held  to  obtain  the  requirements  for  each  analysis.  At  these 
meetings,  the  specific  time  frame,  durations,  and  threshold  levels  were 
specified.  The  type  of  computation  and  form  of  the  output  were  also 
decided.  The  input  parameters  to  the  main  program  were  then  modified  and 
adjustments  to  the  subroutines  were  made  to  provide  the  requested 
computations  and  results.  As  each  record  was  read  from  the  data  file,  it 
passed  through  a  series  of  data  checks.  Records  were  excluded  if  they  were 
not  collected  at  the  Illinois  station,  were  coded  as  bad  data,  or  were 
taken  after  the  last  date  on  the  preceeding  tape.  Castj  were  also  flagged 
if  they  contained  rainfall  durations  that  crossed  seasonal  boundries. 

Once  the  data  had  filtered  through  these  checks,  it  became  part  of  the 
sanitized  data  Let  used  in  the  data  analysis.  The  statistics  routinely 
computed  on  this  file  were: 

*  number  of  years  included  in  the  analysis 

*  number  of  consecutive  minutes  that  pre-determined  time 
intervals  or  threshold  levels  were  exceeded 

*  maximum  duration  rate 

*  average  frequency  count  by  season,  month,  or  year 

*  converting  scores  to  logarithmic  scale 

*  plotting  the  defined  profile 

*  writing  the  results  to  a  separate  file  for  future 
analysis 

The  results  of  each  analysis  were  used  to  formulate  the  specifications  of 
the  next  step  in  the  evaluation  process. 
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2.3.3  Output 


At  the  completion  of  each  run,  the  results  were  summarized  in  tabular 
and/or  graphic  form.  The  objective  was  to  depict  the  output  in  a  way  that 
highlighted  the  salient  features.  The  statistics  were  either  routed  to  a 
printer /plot ter  or  stored  on  a  disk  file  for  future  reference. 

A  list  of  plots  and  tables  generated  by  the  software  is  presented  below: 


PLOTS: 


Seasonal: 


*  overlaid  frequency  plots  of  1,  5,  10,  15,  20,  and  30  minute 
durations  for  up  to  1  nn/min  or  2.5  nn/min. 

*  probability  (Poisson)  of  1,  3,  5,  10,  15,  and  20  minute 
durations  up  to  1.0,  or  2.5,  nn/min. 

*  hourly  interval  (0-4,  4-8,  8-12,  12-16,  16-20,  20-24)  overlaid 
on  1,  5,  10,  15,  20,  and  30  minute  durations  for  up  to 

2.0  nn/min. 


Monthly: 


*  overlaid  frequency  plots  of  1,  5,  10,  15,  20,  and  30  minute 
durations  for  up  to  1  nn/minute. 


Annual: 


*  overlaid  frequency  plots  of  1,  5,  10,  15,  20,  and  30  minute 
durations  up  to  1  nn/minute. 


TABLES: 


Seasonal: 


*  frequency  and  average  intensity  of  rainfall  for  duration 
lengths  ranging  from  1  to  1440  minutes. 

*  year  and  longest  duration  for  thresholds  of  .1,  .2,  .3,  .4, 
.5,  .6,  .7,  .8,  .9,  1.0,  1.3,  1.6,  2.0,  and  2.5  mn/minute 

*  frequency  distribution  if  1  minute  rainfall  rates  by  time  of 
day  (in  hours)  for  varying  threshold  levels. 

*  frequency  of  rainfall  rates  by  hourly  intervals  for 

rates  of  .001,  .25,  .50,  .75,  1.0,  1.5,  and  2.0  mn/minute  and 
rainfall  durations  of  1,  5,  10,  15,  20,  30,  and  35  minutes. 
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Monthly: 


*  year  and  longest  duration  for  thresholds  .1,  .2,  .2,  .4,  .5, 

.6,  .7,  .8,  .9,  1.0,  and  2.5  rar/min. 

Seasonal 

Monthly 

Annual: 

*  frequency  distribution  of  1,  5,  10,  15,  20,  25,  30,  and  35  min 
durations  for  varying  thresholds  from  .01  to  10. 

2.4  Data  File  Description 

2.4.1  Magnetic  Tap  Files 

The  Illinois  data  resided  on  two  magnetic  tapes.  These  tapes  were  both 
internally  unlabeled  and  written  in  fixed  block  formats.  The  first  tape, 
TO=KRAIN1,  contained  35.5  months  of  Illinois  rainfall  data  recorded  between 
May  1969  and  April  1972.  The  second  tape,  TTJ=RAINl,  had  92  months  of 
Illinois  data  collected  from  January  1972  to  August  1979.  A  copy  of  each 
file  was  created  on  disk.  Appendix  A  contains  a  more  detailed  description 
of  the  tape  characteristics. 

2.4.2  Record  Format 

Appendix  B  and  C  contains  the  format  of  the  KRAINl  and  RAINl  data  files, 
respectively.  Inspection  of  these  files  will  reveal  that  the  first  seven 
fields  in  both  files  were  identical.  The  first  field  was  coded  to  either 
indicate  the  beginning/end  of  the  file  or  the  integrity  of  the  data  within 
that  record.  The  second  field  identified  the  station  where  the  data  had 
been  recorded  and  the  third  through  seventh  field  showed  the  year,  month, 
day,  hour,  and  minute  that  the  data  was  recorded.  The  remaining  fields  on 
each  record  were  used  to  record  the  rainfall  rates  to  the  nearest  .001  mm. 
The  KRAINl  file  contained  60  1-minute  measures  of  rainfall  on  each  record 
while  the  RAINl  file  consisted  of  15  1-minute  rainfall  recordings. 
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3.0  DATA  ANALYSIS 

3.1  QYfiiaiew 

This  section  discusses  the  evaluation  process  undertaken  in  this  investi¬ 
gation.  It  is  divided  into  two  parts.  The  first  part  reviews  the  descrip¬ 
tive  statistics  collected  on  the  entire  data  set.  The  second  part  reports 
the  annual,  monthly,  and  seasonal  breakdown  of  the  file.  The  rationale  for 
each  analysis  and  a  synopsis  of  the  results  are  presented. 

3.2  Descriptive  Statistics 

The  first  steps  taken  in  this  analysis  were  to  confirm  the  general  features 
of  the  data.  Using  the  supporting  documentation  provided  with  the  tapes 
and  published  information  about  the  data  base,  a  series  of  descriptive 
statistics  were  developed  to  verify  the  data.  From  these  analyses,  it  was 
determined  that  the  file  contained  5,256,000  1-minute  rainfall  rates  that 
spanned  over  a  ten  year  period.  The  beginning  and  end  dates  were  May  1, 
1969  and  August  31,  1979,  respectively.  Unfortunately,  May  fell  in  the 
middle  of  spring  and  inclusion  of  the  May  1969  data  would  have  artifically 
distort  the  seasonal  statistics.  To  resolve  this  problem,  the  May  1969 
month  was  removed  from  any  analyses.  The  station  identification  number  was 
then  examined  to  insure  that  all  recordings  were  collected  at  the  Illinois 
location.  A  missing  case  report  was  also  performed  to  determine  the  amount 
of  lost  information.  This  report  revealed  that  approximately  109  records 
or  0.5%  of  the  observations  were  unusable. 

The  next  step  was  to  examine  the  frequency  distribution  of  rainfall  rates 
between  .001  and  10  mr/min.  The  results  of  this  analysis  are  presented  in 
Appendix  D.  Inspection  of  this  table  will  show  that  the  increment  between 
various  rainfall  rates  differed.  For  rain  frequencies  between  .001  and 
1.0,  the  gradient  is  .01  whereas  this  value  increases  to  .1,  .5,  or  1.0 
with  higher  intensities.  Since  the  majority  of  rates  occurred  between  .001 
and  1.0,  this  strategy  was  used  to  discriminate  between  sampling  points. 
Rainfall  occurrence  was  most  frequent  between  .001  and  0.5  mr/min  (169108 
observations  in  10.25  years) .  Rains  intensities  drop  off  by  more  than  one 
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order  of  magnitude  in  the  bin  interval  .05  to  .10  im/min  (24069 
observations)  and  interval  .10  to  .15  ira/min  contains  7197  observations, 
after  which  there  is  a  gradual  linear  decline  in  observations. 


3.3  Annual.  Monthly,  and  Seasonal  Analyses 

The  results  of  the  previous  analyses  identified  the  global  characteristics 
of  the  data.  To  obtain  a  better  understanding  to  the  trends  that  existed, 
the  focus  of  analyses  turned  toward  annual,  monthly,  and  seasonal  profiles. 
Rain  intensities  between  0.1  and  1.0  nny'min  (or  2.5  nr/ndn)  with  durations 
of  1,  5,  10,  15,  20,  and  30  minutes  were  evaluated.  Initially  the  inten¬ 
sity  plots  were  broken  down  by  annual,  monthly,  or  seasonal  divisions. 
There  is  little  empirical  information  available  to  predict  the  yearly  or 
monthly  trend  of  rainfall  although  previous  research  has  found  the  greatest 
rates  of  precipitation  in  summer  and  the  least  in  the  winter.  The  outcome 
of  the  seasonal  analysis  supported  these  earlier  findings  (Figures  3.3.1  - 
3.3.4).  There  were  higher  frequency  and  intensity  of  rainfall  for  all 
duration  lengths  in  the  summer  than  in  the  winter.  Spring  and  fall 
duration  frequencies  were  greater  than  winter  but  less  than  summer.  The 
annual  and  monthly  trends  were  not  as  pronounced. 

In  an  effort  to  further  identify  the  nature  of  annual,  monthly,  and 
seasonal  rain,  the  duration  of  selected  rates  were  compared  using  an 
overlaid  plotting  structure.  The  seasonal  results  were  plotted  and 
displayed  in  Figures  3.3.1  through  3.3.4.  Once  again,  the  most 
differentiation  was  found  between  seasons.  At  this  stage  of  the  analysis, 
modifications  were  made  to  the  plot  routines.  Options  were  added  to  the 
program  parameters  to  define  the  length  of  the  X-axis,  convert  abscissa  and 
ordinate  axes  into  logarithmic  scale,  or  choose  the  form  of  the  plots 
(e.g., black  pen,  red  pen). 

ftiset  and  cessation  of  rainfall  is  variable.  It  may  begin  in  one  month  or 
season  and  end  in  a  subsequent  one.  This  phenomenon  can  obscure  the 
boundaries  between  January  and  February  or  Winter  and  Spring.  To 
investigate  this  potential  confound,  the  number  of  rainfall  durations  that 
crossed  montidy  or  seasonal  boundaries  were  examined.  The  frequency  of 
crossovers  were  four  in  forty  one  seasons  or  approximately  10%. 
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SUMMER  PLOT  FOR  SIX  DURATION  TIMES 


RATE  (MM/MI N) 
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The  crossovers  were  handled  by  placing  the  observation  in  the  month  or 
season  where  the  mai^-ity  of  rainfall  occurred. 

The  scope  of  the  invest igation  was  then  directed  toward  the  analysis  of 
rainfall  durations  for  monthly  and  seasonal  categories.  To  determine  more 
specifically  if  there  was  any  clustering  of  long  duration  lengths  by 
intensity  levels,  time  frame  (seasonal  or  monthly),  or  year;  tables  were 
constructed  that  listed  the  longest  duration  length  the  year  that  it 
occured,  and  the  threshold  value  that  it  exceeded.  The  seasonal  table  is 
presented  in  Appendix  E.  The  data  indicate  that  by  season,  the  longest 
duration  of  rates  greater  than  0.1  mr/min  is  in  the  fall,  followed  by 
spring,  sunnier  and  winter. 

The  next  step  was  to  find  the  probability  of  rainfall  occurring  at 
different  intensities  and  durations.  Since  the  frequency  of  rainfall  was 
rare  for  intensities  over  2.5  mny'min  and  durations  greater  than  30  minutes, 
a  Poisson  distribution  was  used.  The  probability  of  obtaining  1,  3,  5,  10, 
10,  or  20  occurences  of  rainfall  that  exceeded  intensities  of  .1  through 
2.5  mrr/min  were  plotted  for  5,  10,  15,  20,  and  30  minute  durations  and  each 
season.  Examples  of  some  plots  are  displayed  in  Figures  3.3.5  and  3.3.6. 

Although  duration  levels  between  1  and  30  minutes  were  mostly  examined  in 
this  study,  continuous  rainfall  did  reach  limits  of  1440  minutes.  Using 
the  intensity  measures,  which  ranged  between  .001  and  10.5,  it  was  possible 
to  assess  the  relationship  between  rainfall  intensities  and  durations  for 
each  season.  There  may  indeed  exist  an  asymptotic  duration  where  rainfall 
in  most  intense.  The  findings  of  this  analysis  showed  that  no  significant 
patterns  were  apparent. 

The  final  breakdown  of  rainfall  was  done  by  time  of  the  day.  TVo  such 
analyses  were  conducted.  The  first  presented  rainfall  intensities  of  .001, 
.25,  .50,  .75,  1.00,  1.50,  and  2.0  im/min  for  each  individual  hour  from 
0000  hours  to  2400  hours.  The  second  analysis  aggregated  these  units  into 
6  4-hour  intervals.  Figures  3.3.7  -  3.3.9  exemplifies  some  seasonal  hourly 
interval  plots. 
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PROBABILITY 


SUMMER  PLOT  -  AT  LEAST  5  OCC. 


o 


Figure  3.3.6 
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4.0  CONCLUSIONS 


A  model  for  estimating  1-minute  rainfall  rates  has  been  developed  using 
stepwise  multiple  regression  analysis.  It  is  made  up  of  six  regression 
equations  to  estimate  rates  that  are  equalled  or  exceeded  0.01,  0.05,  0.10, 
0.50,  1.0,  and  2.0  percent  of  the  time  during  a  month  at  a  given  location. 
The  model  is  based  on  122  months  of  data  gathered  from  12  stations  (Tattel- 
man  and  Scharr  (1983) ) . 

Tahles  and  graphs  presented  in  this  paper  were  used  as  part  of  a 
presentation  given  at  the  4th  Conference  on  Applied  Climatology,  Scotsdale, 
Arizona,  May  21-  24,  1984. 

Hie  Urbana  data,  extensively  examined  here,  is  the  largest  continuous, 
accurate  and  complete  gathering  of  1-minute  instantaneous  rainfall  rates 
ever  recorded.  Hie  analyses  that  have  been  done  on  this  data  base  include 
detailed  descriptive  statistical  information  that  will  be  used  as  a  basis 
for  comparison  with  other  rainfall  data  sources.  The  software  programs 
described  in  this  report  have  been  written  with  flexibility  to  allcw  data 
from  other  sources  and  other  formats  to  be  processed  similarly. 

Further  effort  is  likely  to  include  a  modeling  scheme  where  the  estimate  of 
frequency  of  occurrence  of  rain  durations  for  varying  rain  intensities  is 
done.  Hie  modeling  equation  will  most  likely  be  in  non-linear  exponential 
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APPENDIX  A 


Rainfall  Tape  Characteristics 


Tape  Name 

KRAIN1 

RAINl 

Label 

Unlabeled 

Unlabeled 

Density  (BPI) 

800 

1600 

Record  Length: 

Logical 

253 

80 

Physical 

253 

800 

Record  Format 

Fixed  block 

Fixed  block 

Track  type 

7 

9 

Character  code 

BCD 

ASCI 

No.  of  files 

6 

24 

APPENDIX  B 


Tapel  =  RAINl  (File  24) 
File  Description 


Column  Ha*. 

Format 

Variable  Uame 

identification 

l 

Fl.0 

Code 

0  =  regular  data  record 

1  =  zero  data  record 

2  =  missing  data  record 

3  =  start  record 

4  =  end  record 

2-3 

F2.0 

Station  No. 

TVo  digit  station  ID 
(Urbana,  Ill  =  11) 

4-5 

F2.0 

Year 

Year  of  recording 

72  -  79 

6-7 

F2.0 

Month 

Month  of  recording 

01  -  12 

8-9 

F2.0 

Day 

Day  of  recording 

01  -  31 

10-11 

F2.0 

Hour 

Hour  of  recording 

00  -  23 

12-13 

F  2.0 

Minute 

Minute  of  recording 

00,  15,  30  or  45 

14-73 

15F4.3 

1-minute  rates 

15  1-min  rates  (im/min) 
Rounded  to  nearest  .001. 
Rates  less  than  or  equal 

to  .0042  set  to  0.0. 

appendix  c 


Tape2  =  KRAINl  (file  1) 
File  Description 


Column  No.  Format  Variable  Name 
1  F1.0  Code 

2-3  F2.0  Station  No. 

4-5  F2.0  Year 

6-7  F2.0  Month 

8-9  F2.0  Day 

10-11  F2.0  Hour 


identification 

0  =  regular  data  record 

1  *  zero  data  record 

2  =  missing  data  record 

3  =  start  record 

4  =  end  record 

TVo  digit  station  ID 
(Urbana,  Ill  =  11) 

Year  of  recording 
69  -  72 

Month  of  recording 
01  -  12 

Day  of  recording 
01  -  31 

Hour  of  recording 
00  -  23 


12-13  F2.0  Minute  Minute  of  hour 

00 

14-253  60F4.3  1-min  rates  60  1-min  rates  (mn/min) 

rounded  to  nearest  .001. 
Rates  less  than  or  equal 
to  .0042  set  to  0. 
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APPENDIX  E 


Longest  Seasonal  Duration  Lengths 


Threshold 

(im/min) 

Winter 
length  year 

Spring 

length  year 

£UZHD£I 

length  year 

£all 

length 

year 

0.1 

48 

1971 

75 

1976 

69 

1970 

91 

1969 

0.2 

18 

1975 

45 

1976 

36 

1970 

38 

1969 

0.3 

13 

1975 

39 

1976 

33 

1975 

21 

1978 

0.4 

12 

1975 

37 

1976 

32 

1975 

21 

1978 

0.5 

12 

1975 

32 

1976 

26 

1975 

20 

1978 

0.6 

10 

1975 

15 

1974 

26 

1975 

20 

1978 

0.7 

4 

1975 

15 

1974 

25 

1975 

19 

1978 

0.8 

3 

1975 

14 

1974 

23 

1975 

16 

1978 

0.9 

2 

1976 

14 

1974 

22 

1975 

15 

1977 

1.0 

2 

1976 

13 

1974 

20 

1975 

15 

1977 

1.3 

1 

1975 

12 

1974 

14 

1975 

9 

1977 

1.6 

- 

- 

10 

1974 

10 

1977 

6 

1976 

2.0 

- 

- 

4 

1974 

5 

1977 

3 

1976 

2.5 

— 

— 

4 

1974 

3 

1972 

2 

1976 

82 


i. 


GRAVITY  ANALYSIS 


1.0  INTRODUCTION 
1.1  Backsiound 

The  Air  Force  Geophysics  Laboratory  (AK3L)  is  developing  a  program  to 
verify  gravity  models  by  directly  measuring  gravity  using  high-altitude 
balloons.  Gravity  is  measured  using  a  very  accurate  vibrating  string 
accelerometer  (VSA)  gravimeter  as  part  of  a  balloon  gondola  package. 

Although  the  gravimeter  package  senses  the  acceleration  due  to  gravity,  it 
also  senses  accelerations  due  to  motion  of  the  gondola  package.  These 
accelerations  consist  of  translational  accelerations  of  the  gondola  package 
as  well  as  accelerations  die  to  rotational  motion  of  the  gondola  package, 
i.e. ,  pendulum  effects  and  rotations  of  the  gondola  package  about  its  own 
axis. 

The  pendular  motion  of  the  gondola  package  also  perturbs  the  gravity 
measurements  in  another  way  in  that  the  deflection  of  the  gondola  package 
from  the  local  vertical  (the  direction  of  the  gravity  vectc  )  causes  errors 
in  the  estimate  of  gravity  at  that  point,  since  the  gravimeter  axis  is  no 
longer  aligned  with  the  local  vertical. 

Another  important  consideration  in  using  this  technique  to  estimate  gravity 
is  that  accurate  knowledge  of  the  gravimeter's  position  is  necessary  to 
compare  measured  results  with  gravity  models  at  the  gravimeter  location. 

The  long-term  goal  of  this  project  is  to  determine  gravity  to  1  mgal 
accuracy.  This  "error  budget"  can  be  simply  described  in  terms  of  the 
acceleration,  velocity,  and  position  of  the  gondola  package  and  the 
deflection  from  the  local  vertical  of  the  gondola  package.  Acceleration  of 
the  gondola  package  along  the  local  vertical  of  10~5  n/sec2,  or  deflections 
of  0.08  degrees  from  the  local  vertical  of  the  gondola  package  would  cause 
gravimeter  errors  of  1  mgal.  Also,  a  correction  to  the  gravimeter 
measurement  rrust  be  made  that  depends  on  the  east -component  of  the  gondola 
package's  velocity.  An  error  of  0.08  n/sec  in  the  estimate  of  this 
velocity  would  cause  an  error  in  the  gravity  estimate  of  1  mgal.  Also, 


since  gravity  varies  in  the  vertical  direction  roucjily  as  1  mgal  per  3 
meters,  errors  in  the  estimate  of  the  gondola  package  of  3  meters  would 
also  cause  gravity  estinate  error  of  1  mgal. 

Given  that  the  actual  acceleration,  deflection,  horizontal  velocity  and 
position  errors  are  orders  of  magnitude  larger  than  the  limits  shewn,  using 
unprocessed  gravimeter  data  as  the  gravity  estimate  would  be  completely 
unsatisfactory.  The  solution  proposed  by  APGL  involves  the  accurate 
determination  of  the  motion  of  the  gondola  package  (position,  velocity, 
acceleration,  and  rotation)  and  use  of  this  information  to  correct  the 
gravimeter's  measurement.  The  combination  of  an  accurate  tracking  system 
and  an  inertial  motion  sensing  package  (IMSP)  is  required  in  order  to 
estimate  gondola  position,  velocity,  and  acceleration  accurately. 

In  October  1983,  a  launch  of  a  high-altitude  balloon  was  conducted  by  APGL 
personnel  at  Holloman  AFB  in  New  Mexico.  The  balloon  gondola  package 
included  a  VSA  gravimeter  and  an  inertial  motion  sensing  package  which 
included  3  orthogonal  accelerometers,  3  orthogonal  rate  gyros,  and  3 
orthogonal  magnetometers.  Gondola  position  tracking  was  done  by  three 
ground-based  digital  radars  located  at  White  Sands  Missile  Range  (WSMR) . 

The  objectives  of  this  flight  were  to  demonstrate  the  feasibility  of 
measuring  gravity  using  a  high-altitude  balloon  and  to  evaluate  the 
performance  of  the  digital  radar  in  tracking  the  balloon.  The  objective  of 
this  flight  was  not  to  estimate  gravity  to  1  mgal  accuracy,  although  it  was 
hoped  that  the  data  could  be  processed  in  some  way  to  improve  the  raw 
gravimeter  output. 

1.2  Overview 

BEDFORD  RESEARCH'S  participation  in  the  post-flight  data  analysis  included 
the  creation  of  an  investigator's  data  tape  and  the  support  of  AEGL 
personnel  in  their  analysis  of  gravimeter,  tracking,  and  IMSP  data.  The 
investigator's  data  tape  is  a  data  base  consisting  of  evenly-sampled  data 
from  the  VSA  gravimeter,  IMSP,  digital  radar,  and  other  sources. 


84 


The  following  3  sections  of  this  report  will  describe  3  separate  technical 
projects  involving  the  processing  of  the  balloon  gravity  data.  Section  2.0 
describes  an  algorithm  used  to  predict  missing  data  found  in  the  radar 
tracking  channels.  Section  3.0  describes  a  computer  program  which  allows 
users  to  read  selected  time  segments  and  channels  of  the  investigator's 
data  tape  interactively  and  to  place  these  data  into  separate  I*S*P  data 
files  for  subsequent  processing  and  analysis.  Finally,  Section  4.0, 
describes  sane  analysis  of  the  radar  tracking  data  which  was  intended  to 
get  a  better  estimate  of  the  balloon  trajectory.  This  analysis  involved 
the  use  of  a  Kalman  filter  to  incorporate  radar  tracking  data  and  VSA  data 
in  order  to  smooth  the  vertical  position  estimate  of  the  gondola  location. 


2.0  PREDICTION  FILTER  TO  FILL  IN  GAPS  OF  RADAR  TRACKING  DATA 
2.1  Introduction 

One  of  the  problems  encountered  during  the  post-flight  data  analysis  was 
the  existence  of  numerous  data  gaps  in  the  WSMR  radar  tracking  data.  The 
presence  of  gaps  in  the  data  can  cause  severe  problems  for  tracking 
algorithms,  and  may  have  caused  problems  for  DMA  in  deriving  gravity 
estimate  values  from  tracking  data.  For  these  reasons  an  algorithm  was 
developed  to  predict  the  missing  data  to  preserve  its  continuity. 

Sixty-or»e  gaps  were  identified  which  ranged  from  2  to  30  samples  in  length. 
In  addition,  three  large  gaps  were  identified  which  were  259,  19S,  and  839 
samples  long.  Because  of  their  length,  these  gaps  were  not  filled  since  the 
performance  of  any  algorithm  would  be  highly  suspect  over  such  long 
periods. 

This  section  presents  the  algorithm  used  to  predict  the  missing  data.  This 
algorithm  is  based  on  the  theory  of  linear  prediction  [1] .  Following  a 
discussion  of  the  algorithm,  a  brief  discussion  of  its  implementation  and 
some  graphical  results  are  presented. 
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2.2  Linear_Prediction  Algorithm 

The  approach  taken  to  fill  in  the  gape  is  to  use  the  technique  of  linear 
prediction.  This  technique  assumes  an  auto-regressive  (AR)  time  series 
model  for  the  data,  and  uses  this  model  to  predict  the  missing  data  based 
on  the  existing  data  surrounding  the  gap.  The  signal  x(n)  is  assumed  to  be 
a  linear  combination  of  past  values  and  some  input  u(n): 

P 

x(n)  =  ^  \x(n-k)  +  u(n)  (1) 

k=l 


The  input  u(n)  is  assumed  to  be  a  white  noise  term  with  variance  <*  2.  The 
model  order  P  must  be  estimated  and  the  AR  coefficients  {a^}  must  be 
determined.  Experiments  with  the  tracking  data  showed  that  P=8  is  a  good 
model  order  to  use. 

The  sample  values  x(nQ)  f  x(no+l) ,. . .  ,x(ni)  are  assumed  to  be  missing. 
Thus,  the  gap  is  L-n^+i  points  wide. 

The  data  immediately  preceding  the  gap  are  used  to  estimate  the  forward 
prediction  coefficients  {a^j  via  ^e  Burg  method  [2] .  This  method  is  part 
of  the  basic  i*S*P  package  under  the  function  name  LNP.  500  points  are 
used  to  form  the  estimates.  These  prediction  coefficients  are  used  in  the 
following  formula  to  yield  the  forward  prediction  data: 


xf  (n) 


^  akxf(n-k),  n=no,no+l,...,ni 


(2) 


where  f0r 

Similarly,  the  data  iirmediately  following  the  gap  are  used  to  estimate  the 
backward  prediction  coefficients  Again,  the  Burg  method  is  used  and 

500  points  are  used  to  form  the  estimate.  The  backward  prediction 
coefficients  are  used  in  the  following  formula  to  form  the  backward 
prediction  data: 


p 


bkXbto+k)*  n=ni,ni-l,...,no 


(3) 


where  3^)(n}=x(n)  for  n>ni. 

The  final  prediction  estinate  x(n)  is  formed  by  adding  a  weighted  version 
of  the  forward  and  backward  estimates.  The  weights  are  proportional  to  the 
error  in  the  predictors  where  the  error  will  increase  as  n  is  farther  away 

from  nQ  for  the  forward  predictor  and  n^  for  the  backward  predictor.  The 
final  prediction  equation  is  then  given  by: 

x(no+m)  =  (n^r5-)  Xf(no+m)  +  ^^Xb(no+m) ,  (4) 

m  ■  0,  1 ,  . .  * ,  L- 1 


2.3  Implementation  Issues 

If  the  data  have  a  DC  offset  anchor  linear  trend,  they  can  cause  errors  in 
estimating  the  prediction  coefficients  and  will  adversely  affect  the 
prediction  results.  Therefore,  the  DC  offset  and  trend  re  first  removed 
from  the  data  preceding  the  gap  and  from  the  data  following  the  gap  before 
estimating  the  prediction  coefficients.  The  forward  and  backward 
prediction  coefficients  are  then  computed  from  the  detrended  data.  The 
detrended  forward  and  backward  predictions  are  then  computed  using  (2)  and 
(3).  Ihe  DC  offset  and  trend  are  then  added  back  to  the  predictions  and  the 
weighting  in  (4)  is  finally  applied. 

One  problem  that  was  encountered  during  the  graphical  inspection  of  the 
gaps  was  that,  very  often,  the  data  immediately  surrounding  the  gaps  seemed 
to  be  suspect.  This  is  shown  in  Figs.  2. 1-2. 3  where  the  gap  is  indicated 
ty  the  two  vertical  lines.  Note  that  Figs.  2.2  and  2.3  shew  two  gaps  in 
the  time  span  shewn.  These  figures  shew  that  the  data  just  before  the  gap 
and  just  after  the  gap  are  noisier  than  data  further  away  from  the  gap. 
For  this  reason  it  was  decided  to  extend  the  size  of  the  gap  (or  to  predict 
note  points)  for  gaps  that  shewed  tiers  tendency.  Also  displayed  in  Figs. 
2. 1-2.3  ate  the  pt  ei  lift  ions  ioi  the:*-  gaps.  Hie:*'  ai  v  the  junooth  turves 
that  joss  over  the  gaps.  In  Figs.  2.2  and  2.3  wliere  two  gaps  are  so 
closely  located,  the  two  gaps  ore  treated  as  a  single  gap. 
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3.0  PROGRAM  TO  READ  INVESTIGATORS'  DATA  TAPE 

3.1  Introduction 

The  investigators'  data  tape  (IOT)  is  a  unified  data  base  containing  48 
channels  of  data  plus  time  code.  The  data  were  put  on  the  tape  at  a  sample 
rate  of  10  Hz. 

This  section  briefly  describes  a  FORTRAN  program  written  to  allow  users  to 
interactively  read  selected  time  segments  and  channels  of  this  file  and 
place  the  data  into  separate  I*S*P  data  files  for  subsequent  processing  and 
analysis. 

3.2  ProgranLDescription 

A  FORTRAN  program  was  written  to  allow  users  to  interactively  access  the 
IDT  and  place  particular  channels  and  time  segments  into  I*S*P  files  for 
further  analysis. 

The  program  will  prompt  the  user  for  the  necessary  input  and  output 
information  and  then  open  the  IDT  (which  is  actually  a  disk  file  on  the 
VAX) ,  read  the  appropriate  channels,  and  write  the  specified  I*S*P  files. 

The  program  begins  with  the  following  prompts: 

Enter  start  time  (Zulu)  below  : 

Enter  hours  : 

Enter  minutes  : 

Enter  seconds  : 

The  user  should  type  in  the  time  at  which  to  begin  reading  the  IOT  as 
integers  by  first  specifying  the  time  in  hours,  then  minutes,  then  seconds. 
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The  user  will  then  be  prompted  for  the  amount  of  time  he  wishes  to  read  in 
the  same  way  he  was  prompted  for  the  start  time: 

Enter  amount  of  time  to  read  below: 

Enter  hours  : 

Enter  minutes  : 

Enter  seconds  : 

Next,  the  following  prompt  appears: 

Enter  decimation/averaging  factor  : 

The  user  should  enter  1  if  he  wishes  to  preserve  the  original  10  Hz 
sampling  rate,  or  an  integer  greater  than  1  if  he  wishes  to  decimate  or 
average  the  data.  Caution  should  be  taken  when  choosing  a  factor  greater 
than  1  since  aliasing  will  occur.  The  proper  way  to  perform  sampling  rate 
decimation  is  to  choose  a  decimation  factor  of  1  and  then  use  the  I*S*P 
program  SRC  to  decimate  or  interpolate  the  files  to  the  desired  output 
sampling  rate  using  the  proper  decimation/ interpolation  techniques  [3]. 

If  a  decimation  factor  other  than  1  is  chosen,  the  following  prompt 
appears: 

Enter  A  for  average  or  D  for  decimation  : 

If  A  is  selected,  n  samples  will  be  averaged  to  form  an  output  value  and 
the  output  samples  will  be  separated  by  0.1n  sec,  where  n  is  the 
decimation/averaging  factor.  If  D  is  selected,  every  nth  sample  will  be 
output  and  again,  the  output  samples  will  be  separated  by  0.1n  sec. 

The  program  will  then  inform  the  user  how  many  records  will  be  written,  the 
record  size,  the  amount  of  time  to  be  read,  and  will  ask  if  that  is 
acceptable  by  typing  the  following  prompt: 

Is  that  OK  (Y  or  N) ?  : 
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The  user  should  respond  with  a  Y  for  yes  or  an  N  for  no.  If  Y  is  typed, 
the  program  will  proceed  and  if  n  is  typed,  the  program  will  go  back  to  the 
beginning  and  prompt  for  the  starting  time  at  which  point  the  user  can  type 
in  new  values. 

The  next  prompt  is: 

Enter  number  of  channels  to  write  : 

The  user  can  then  respond  with  an  integer  from  1-48. 

The  next  prompt  is: 

Enter  channel  number  (0  for  menu)  : 

The  program  will  type  this  prompt  from  1  to  47  times,  depending  on  the 
number  of  channels  entered  in  the  previous  prompt.  If  48  is  entered,  all 
channels  will  be  written  so  the  user  is  not  prompted  for  channel  numbers. 
If  a  0  is  typed,  a  menu  of  the  48  channels  is  typed  on  the  screen. 

The  last  prompt  is: 

Enter  first  ISP  file  number  for  output  : 

The  user  should  enter  an  integer  from  0-999.  In  the  current  version,  if 
more  than  one  channel  is  written,  an  I*S*P  file  is  written  for  each  channel 
and  the  numbering  is  consecutive  starting  at  the  number  typed  in  for  the 
first  file.  In  future  versions  this  may  be  changed  to  take  advantage  of 
I*S*P's  multichannel  file  capability. 

The  following  is  a  sample  output  from  the  program. 

Enter  start  time  (Zulu)  below  : 

Enter  hours  :  17 
Enter  minutes  :  52 
Enter  seconds  :  0 

Enter  amount  of  time  to  read  below: 


CV! 
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Entei  hours  :  0 
Enter  minutes  :  10 
&iter  seconds  :  0 

Enter  decimation/averaging  factor  :  1 
Records  to  be  written t  3 
Record  size:  2048 

Amount  of  time  to  be  read:  0  hr  10  min  14  sec 
Is  that  OK  (Y  or  N) ?  :  Y 
Enter  number  of  channels  to  write  :  2 
Enter  channel  number  (0  for  menu)  ;  29 
Enter  channel  number  (0  for  menu)  :  38 
alter  first  ISP  file  lumber  for  output  :  900 
1 


Number  of  records  written* 
FORTRAN  STOP 


4.0  KALMAN  FILTER  TO  SMDOTH  TRACKING  ESTIMATES 


4.1  introduction 


Inspection  of  WSMR  tracking  data  and  errors  raised  some  questions  as  to  the 
validity  of  the  error  estimates.  WSMR  reported  tracking  errors  on  the 
order  of  3  m  for  iruch  of  the  flight,  but  visual  inspection  of  these  data 
along  with  an  understanding  of  the  basic  physics  of  the  balloon  gondola 
package  imply  that  the  dynamics  seen  in  the  data  are  unreasonable  and 
unrealistic.  It  was  therefore  postulated  that  the  errors  as  reported  by 
WSMR  are  overly  optimistic. 

To  try  to  get  a  better  estimate  of  the  balloon  trajectory,  a  Kalman  filter 
was  designed  to  incorporate  both  tracking  data  and  accelerometer  data  from 
the  gondola  package.  Section  4.2  is  devoted  to  a  description  of  the  model 
and  assumptions  used  to  yield  the  estimate  of  the  balloon  trajectory,  and 
Section  4.3  shows  the  results  of  the  Kalman  filter  output. 


4-2  Kedman_FiItfii_Design 


For  simplicity  only  the  vertical  direction  was  considered.  Hie  radar 
tracking  data  had  been  converted  to  geodetic  coordinates  so  the  geodetic 
height  as  determined  by  WSMR  radar  was  used  for  the  analysis.  It  was  also 
assumed  that,  at  float,  the  balloon  pendulous  motion  was  small  so  that  the 
z-accelercmeter  and  VSA  were  reasonably  aligned  with  the  local  vertical,  or 
geodetic  height  direction.  Hie  gondola  geodetic  height,  h,  velocity, 
dh/dt,  and  acceleration,  d^h/dt2  are  modelled  in  a  three  state  Kalman 
filter  model  as: 


x(k+l)  =  x(k)  +  Gu(k) 


x(k+l)  = 


1 

0 

0 


A 

1 

0 


.5A2 

A 

a 


x(k) 


u(k) 


where  x(k)=(x^(jtj  x2(k)  x3(k))T'  xl(k)=h,  X2(k)=dh/dt,  x3(k)  =d^h/dt2,  and 
=0.1  sec  is  the  sampling  interval.  For  the  data  sequent  used  (a  ten  minute 
period  starting  at  1752  Zulu),  the  following  parameters  are  valid: 


01  =  0.994 

q  =  E[u2(k)]  =  1.42  x  10"8. 

The  observations  used  are  the  geodetic  height  as  determined  by  WSMR  and  the 
VSA  output.  Hie  geodetic  hei^it  was  first  pre-processed  by  subtracting  the 
value  of  the  first  point  so  that  the  observations  would  start  at  zero.  Hie 
VSA  data  was  pre-processed  by  removing  the  mean  value  (due  mainly  to 
gravity),  and  by  converting  to  units  of  it/ sec2. 
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In  the  Kalman  filter  formalism  the  observations  are  written  as 


z(k)  =  1  0  0  x(k)  +  v(k) 

.  0  0  1  . 


where  z(k)  and  v(k)  are  (2  X  1)  dimension  vectors  representing  the 
observation  vector  and  the  observation  noise  vector  respectively.  The 
observation  noise  has  a  covariance  matrix  defined  by 

R  =  E[v(k) v(k)T] 

=  r  100  0 
_  0  2.5E-7  . 

Note  that  the  values  used  here  imply  a  standard  deviation  of  10  m  for  the 
tracking  data  and  5E-4  n/sec^  for  the  VSA  data.  The  10  m  standard 
deviation  for  the  tracking  data  is  roughly  10  times  larger  than  the  values 
given  by  WSMR. 


4.3  Results 

Using  the  model  developed  and  the  parameters  shown  in  Section  4.2,  the 
Kalman  filter  was  implemented  to  obtain  estimates  of  the  state  vector  x(k) 
for  the  time  period  shown  in  Fig.  4.1.  The  result  for  the  geodetic  height 

^=xl(k)  is  shown  in  Fig.  4.2.  The  ragged  curve  is  the  original 
observations  as  seen  in  the  top  curve  in  Fig.  4.1  and  the  smooth  curve  is 

the  Kalman  filter  output  estimate  for  the  geodetic  height  h=x1(k)> 

These  results  offer  more  evidence  that  the  errors  given  by  WSMR  are  overly 
opimistic.  Noting  the  excursions  in  the  origiral  observations,  it  is 
clear  that  the  errors  are  much  greater  than  the  3  m  reported  by  WSMR. 
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Figure4,l  -  Observations  used  as  inputs  to  the  Kalman  filter.  Top  curve  is  the 
geodetic  height  in  meters  with  the  initial  value  subtracted  off.  Initial  value 
is  30218.2  m.  Second  curve  is  VSA  data  in  units  of  m/sec^  with  the  mean  value 
subtracted.  Time  units  are  seconds  with  a  start  time  of  1752  Zulu. 
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THEORY  AND  MODEL  DC  OF  INFRARED  EMISSION 


1.0  INTRODUCTION 

The  work  consists  of  generating  and  maintaining  computer  codes  and  other 
software  involved  in  the  generation  of  predictive  codes  of  infrared 
emission  from  quiescent  and  disturbed  atmospheres. 

2.0  DATA  ANALYSIS 
2.1  Data  Modeling 

The  first  project  undertaken  was  the  generation  of  code  to  fit  a  synthetic 
spectrum  to  the  nitric  oxide  (NO)  molecular  fundamental  vibration-rotation 
bands  near  5.3  m  m  wavelength  (1890  cnT1  frequency)  by  using  nonlinear 
minimization  techniques.  The  parameters  to  be  fit  were  the  number  densi¬ 
ties  of  molecules  in  the  upper  vibrational  state  (v1)  of  the  transition 
that  is,  nyI  (molecules/cm^)  and  temperature  (T) .  The  fitting  problem  is 
nonlinear  in  temperature  T,  but  linear  in  the  vibrational  populations  ny, _ 

The  quantity  to  be  minimized  (I)  was  the  mean-square  error  between  the 
measured  spectral  radiance  Iexp(  v  )  (photons/cm2sec  sr  cnT1)  (sr  =  star 
radians)  where  is  the  frequency  in  wavenumber  (cm-1) ,  and  the  calculated 
synthetic  spectral  radiance  I  ( *  rf  J f  T) .  $  is  given  by 

/»'  max 

d  ll(  „  ,{nv.|,T)  -  IexpHl2 
*’  min 

The  spectrum  is  sampled  at  0.25  an-1  intervals. 

The  form  of  the  synthetic  spectrum  I  is  as  follows: 


Iu!  =  nyI  g  e  I-C2(E'-Ev.i/21)A1 

4  QV(T) 

where 

A  *=  the  Einstein  A-  coefficient  in  (sec-1) 

V  *  the  population  density  (molecules/an1) 


BE 


i  ryrvTV  '"jrrjrm  ir*ir*v*  VTWW7 


* 


gj '  =  degeneracy  of  the  state  j*  =  2 j '  +1 

c2  -  second  radiation  constant  =  1.4388  k/cnf^ 

E'  =  energy  of  the  upper  level 

^v'1/21  ~  the  energy  of  state  when  =  v' ,  j  =  1/2,  i  =  1 
As  =  an  effective  path  length 

ul  means:  from  the  upper  transition  (v1 j'i')  to  the  lower  transition 
(v"j"i") 

Qyi  (t)  «  partition  sum 


-  n^rot) 

~  Qy.  (T) 


Q(n)(T) 


Qv'fl  (T)  =  1  +  e~°2A/1!  A  =  e2  -  El  *=  120  cm-1 

Qv'(rot)T)  =— —  el3c2Bv'/4T]  By*  =1.6  -1.7  cm-1 

c2Bv' 

To  about  1%,  the  exponential  is  approximately  unity  and  can  be  neglected 
giving 


T 

c2Bv* 
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2.2  Data  Manipulation 

The  next  project  consisted  of  transferring  SPIRE  CVF  data  from  Computer 
Center  {CC)  tapes  to  tapes  with  proper  format  for  use  on  the  APOLLO 
computer  system.  This  involved  reformatting  the  tapes  and  writing  them 
without  labels. 

2.3  Program  Modification 

The  final  project  consisted  of  modifying  code  to  calculate  the  zenith 
auroral  spectral  intensity  for  NO  and  NO*  band  systems.  This  involved 
adding  subroutines  to  the  zenith  auroral  programs  already  running,  which 
calculated  the  spectrum  at  a  given  altitude.  Cnee  completed,  the  results 
of  the  zenith  auroral  spectral  program  are  convolved  with  a  Hamming  window 
line  shape  and  then  plotted. 
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BOREHOLE  LOG  PLOTTING  AND  ANALYSIS 


1.0  OBJECTIVES 

The  objectives  of  this  project  were  first  to  demonstrate  that  a  program 
possessed  by  the  Solid  Earth  Geophysics  Branch  (IWH)  would  reformat 

geophysical  well  data  tapes.  Seoondly,  a  program  to  produce  plots  of  the 
borehole  log  data  using  BEDFORD  RESEARCH'S  graphic  software  was  desired. 
Thirdly,  any  additional  analyses  were  to  be  provided  as  required. 

2.0  ANALYSIS 

2.1  Data  Conversion 

The  program  SCFSCH  converts  geophysical  well  data  written  in  Log 

Information  Standard  (LIS)  format  to  Log  Information  Library  (LIL)  format. 
The  program  was  written  by  Bendix  Field  Engineering  Corporation,  for  the 
CDC  system  and  was  copied  to  the  Cyber  at  the  Air  Force  Geophysics 

Laboratory.  To  run  SCFSCH,  it  was  first  necessary  to  type  in  a  channel 

identification  table  provided  in  the  documentation.  A  program  was  then 
written  to  copy  the  tape  file  to  disk  by  buffering-in/buf lering-out  each 
record.  Program  SCFSCH  can  then  be  run  as  a  batch  job,  attaching  the 
necessary  files,  creating  a  file  in  LIL  format  and  purging  the  original 
file  in  LIS  format. 


The  log  plotting  program  (LOGPLOT)  was  written  for  the  VAX.  The  program 
reads  the  well  log  files  written  in  LIL  format  and  plots  the  data  requested 
by  the  user.  The  program  utilizes  BEDFORD  RESEARCH  graphics  software. 
This  enables  the  user  to  run  LOGFLOT  on  either  a  VT125  or  Tektronix 
terminal. 

There  are  several  options  available  with  LOGELOT.  The  user  selects  the  log 
file  and  chooses  the  number  and  types  of  curves  he  wishes  to  plot.  He  also 
selects  the  depth  ranges,  the  curve  value  ranges,  and  the  scale  factor  for 
each  axis.  Additional  options  include  the  plotting  of  the  center  line  and 
the  selection  of  a  linear  or  logarithmic  curve  axis.  Figures  2.1,  2.2,  and 
2.3  show  3  log  plots  using  the  various  options. 


The  plots  on  the  VAX  are  produced  after  the  log  file  in  LIL  format  has  been 
transferred  from  the  CYBER.  The  plots  can  be  printed  on  a  hard  copy  unit, 
on  8  1/2  x  11  inch  paper,  with  good  resolution,  resulting  in  a 
journal-quality  figure.  However,  for  research  purposes,  the  optimal  depth 
scale  is  2.5  -  5  inches/100  feet.  Since  the  boreholes  range  from  approxi¬ 
mately  2000  to  8000  feet  deep,  this  necessitates  plots  that  are  at  least  50 
inches  long.  It  would  not  be  feasible  to  produce  these  plots  on  the  VAX 
because  there  is  no  hard  copying  unit  available  for  such  lengthy  plots. 
Therefore,  the  Initiator  requested  a  version  of  LOGPLOT  on  the  CYBER.  This 
would  enable  the  use  of  the  Calcomp  plotter,  on  which  there  is  no  limit  to 
the  length  of  the  x  (or  depth)  axis.  The  options  are  similar  to  those 
available  on  the  VAX  program,  except  that  on  the  CYBER  there  is  a  maximum 
of  4  curves  per  plot  since  there  are  only  4  colored  pens. 

2 . 3  Interactive^ Programs 

BEDFORD  RESEARCH  personnel  wrote  a  procedure  file  for  the  CYBER  which 
enables  the  Initiator  to  run  programs  SCFSCH  and  LOGPLOT  easily.  The 
procedure  file  runs  interactively  and  allows  the  user  to  select  one  of 
three  options:  1)  copy  an  LIS  tape  to  a  disk  file,  2)  un  SCFSCH,  and 
3)  produce  a  log  plot.  The  procedure  prompts  the  user  for  the  information 
necessary  to  generate  a  batch  job  for  the  particular  option  selected. 
After  executing  the  procedure  file,  the  batch  job  is  sutmitted. 

IXiring  the  execution  of  program  SCFSCH,  only  part  of  the  information  from 
the  original  LIS  file  is  decoded.  Therefore,  the  Initiator  requested  a 
program  to  edit  the  created  LIL  file  and  incorporate  additional  data  such 
as  the  geographical  coordinates  of  the  borehole,  the  casing  and  logging 
information,  etc.  The  resulting  LILEDIT  program  on  the  VAX  modifies  6 
records  and  adds  5  records  to  the  file  in  a  format  that  meets  LIL 
specifications. 


3.0  SUMMARY 


Research  involving  the  well  log  data  was  also  undertaken  by  BEDFORD 
RESEARCH  personnel.  This  involved  two  separate  projects:  1)  an  analysis  or 
tube  waves  in  the  cased  Zeeck  well  and  2)  a  comparison  of  seismic 
velocities  for  6  wells  in  the  Texas  Panhandle.  Figure  3.1  shws  the  drift 
curves,  the  difference  between  the  integrated  sonic  time  and  seismc  tune 
for  the  wells.  Ihe  preliminary  analysis  for  these  projects  was  completed 
but  the  projects  were  not  brought  to  their  final  completion,  due  to 
complications  and  lack  of  time. 
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Figure  3.1  -  Drift  Curve 
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